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1.  INTRODUCTION 


This  document  is  the  final  report  on  Contract  F19628-80-C-0187.  The  purpose  of 
the  research  reported  herein  is  to  study  the  feasibility  of  determining  ionospheric 
electron  density  profiles  from  measurements  of  the  scattering  of  electromagnetic 
waves  entering  the  ionosphere  from  outside  the  earth's  atmosphere.  This  is  an  inverse 
scattering  problem,  and  the  approach  used  is  to  develop  a  practical  procedure  based  on 
inverse  scattering  theory. 

A  detailed  literature  survey  on  the  application  of  existing  approaches  to  the 
determination  of  the  electron  density  distribution  in  the  ionosphere  led  to  the 
conclusion  that  the  best  approach  was  to  apply  the  Exact  Inverse  Scattering  Theory, 
developed  by  N.  N.  Bojarski,  to  the  problem.  This  literature  survey  is  contained  in 
Volume  II  of  this  report.  Comments  on  the  results  of  the  survey  are  included  in  this 
section.  This  section  also  defines  the  specific  problem  being  addressed  in  more  detail, 
and  examines  the  relationships  among  the  desired  quantity,  measurable  quantities,  and 
the  physical  parameters  useful  from  a  theoretical  standpoint.  Section  2  reviews  the 
derivation  of  the  Exact  Theory.  In  addition,  the  theory  is  cast  into  a  form  including  full 
vector  electromagnetic  fields  and  the  effects  of  a  uniform  magnetic  field.  Section  3 
examines  the  uniqueness  of  the  solution  obtained  using  the  Exact  Theory,  and  extends 
the  author's  earlier  uniqueness  proof  for  the  scalar  case  to  the  full  anisotropic  medium, 
vector  electromagnetic  case.  Prior  to  1980,  the  uniqueness  of  the  solution  to  the 
inverse  scattering  problem  was  complicated  by  the  question  of  nonradiating  sources.  In 
1980,  this  author  presented  a  proof  that  while  nonradiating  sources  are  mathematically 
permissible  within  the  context  of  inverse  scattering  theory,  the  fields  they  produce  are 
physically  inconsistent  with  the  wave  equation  governing  the  problem.  During  the 
course  of  this  research,  two  additional  independent  verifications  that  nonradiating 
sources  are  physically  unrealizable  became  available,  as  reported  in  Section  3. 

Perhaps  the  most  important  and  exciting  results  to  come  out  of  this  research  are 
presented  in  Section  4.  Late  in  1980,  Bojarski  obtained  an  exact  closed-form,  analytic 
solution  for  the  total  fields  (including  fields  in  the  source  region)  in  the  inverse 
scattering  problem.  This  solution  has  the  same  importance  for  the  inverse  problem  as 
the  Kirchhoff  and  Stratton-Chu  equations  have  for  the  direct  scattering  problem. 
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Subsequently,  the  author  was  able  to  obtain  two  results.  First,  the  solution  obtained  by 
Bojarski  was  cast  into  a  form  involving  only  known  quantities  and  Fourier  transforms  of 
the  measured  data.  This  makes  the  new  solution  for  the  fields  practical  for  application 
to  problems  such  as  the  electron  density  profiling  problem.  It  also  renders  the 
computations  associated  with  evaluating  the  new  solution  numerically  as  efficient  as 
the  fast  Fourier  transform,  both  in  speed  and  storage.  Second,  and  even  more 
important,  the  author  was  able  to  obtain  an  exact,  closed-form,  analytic  solution  for 
the  source  term  in  the  inverse  scattering  problem.  This  also  is  expressed  in  terms  of 
only  known  quantities  and  Fourier  transforms  of  the  measured  data.  Perhaps  the  most 
pleasing  aspect  of  these  results  is  the  tremendous  amount  of  new  insight  into  the  basic 
physics  of  inverse  scattering  (and,  indeed,  into  fundamental  field  theory)  which  they 
provide.  As  should  be  the  case  in  a  true  research  effort,  these  results  were  not 
anticipated  before  the  research  was  begun. 

It  is  worth  putting  these  new  results  into  perspective  with  the  previous  status  of 
the  Exact  Theory,  on  which  they  are  based.  Until  now,  it  was  necessary  to  solve  an 
integral  equation  numerically  to  obtain  the  sources,  and  thence  the  refractive  index  and 
the  electron  density  profile.  Although  this  equation  should  have  an  efficient, 
numerically  stable  method  of  solution  it  did  not  have  the  advantages  of  an  analytic 
solution.  The  new  results  provide  analytic  solutions  for  both  the  fields  and  the  sources 
in  the  inverse  scattering  problem. 

Section  5  summarizes  the  results  of  analytical  and  numerical  simulation  work 
performed  to  demonstrate  how  the  electron  density  distribution  can  be  determined 
using  inverse  scattering.  It  was  originally  desired  to  simulate  scattering  by  an 
ionosphere  that  is  inhomogeneous  in  two  or  three  dimensions,  and  apply  the  inverse 
scattering  techniques  to  obtain  a  reconstruction  of  the  simulated  electron  density 
distribution.  Substantial  problems  were  encountered  in  carrying  out  the  simulations, 
because  of  difficulties  in  modeling  the  direct  scattering  problem.  These  results  are 
necessary  because  they  constitute  the  simulated  data  used  as  input  for  the  inverse 
scattering  calculations.  However,  it  was  possible  to  demonstrate  the  stability  of  the 
numerical  solution  for  the  inverse  scattering  approach.  Several  methods  of  computing 
the  required  direct  scattering  data  were  tried,  and  a  number  of  approximations 
heretofore  overlooked  in  the  literature  on  scattering  by  inhomogeneous  media  in 
general  and  on  scattering  by  the  ionosphere  in  particular  were  identified.  These  results 
have  important  implications  for  other  modeling  studies.  A  method  capable  of  providing 
the  required  direct  scattering  data  was  developed,  and  preliminary  tests  indicate  that  it 
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will  provide  the  needed  data  for  the  desired  two-  and  three-dimensional  simulations  in 
future  work.  Finally,  as  means  of  demonstrating  what  should  be  achievable  in  applying 
the  Exact  Theory  to  two-  and  three-dimensional  inhomogeneities,  a  related,  one¬ 
dimensional,  time  domain  simulation  was  carried  out.  The  results  demonstrate  the 
capability  to  reconstruct  realistic  one-dimensional  ionospheric  profiles,  in  the  presence 
of  realistic  experimental  noise.  Furthermore,  the  capability  to  determine  the  profile 
above  the  maximum  in  electron  density  using  only  bottom-side  vertical  incidence 
sounder  data  was  demonstrated.  To  the  author's  knowledge  this  is  the  first  time  such  a 
capability  has  been  demonstrated  even  in  simulation  in  the  presence  of  realistic 
amounts  of  measurement  noise.  The  technique  should  be  directly  applicable  to  the 
analysis  of  data  from  a  large  class  of  existing  ionospheric  sounders. 

Section  6  discusses  the  practical  aspects  of  applying  the  theory  and  results  of  the 
previous  sections  to  the  experimental  determination  of  the  electron  density  distribu¬ 
tion.  In  particular,  it  is  shown  that  while  the  effects  of  the  earth's  magnetic  field  can 
be  included  in  the  analysis,  it  is  not  necessary  to  do  so.  Proper  choice  of  the  frequency 
of  the  electromagnetic  wave  used  to  probe  the  ionosphere  eliminates  this  complication 
and  simplifies  the  required  measurements. 

Section  7  presents  the  conclusions  drawn  from  this  research  and  provides 
recommendations  for  further  work. 

This  paper  has  aspects  of  both  a  report  of  new  results  and  a  review.  This  was 
done  in  a  conscious  effort  to  make  the  result  self-contained.  Readers  familiar  with  the 
author's  earlier  review  papers  on  exact  inverse  scattering  theory  should  feel  free  to  skip 
those  sections  which  appear  familiar. 

1.1  PROBLEM  DEFINITION  AND  THE  CONCEPT  OF  INVERSE  SCATTERING 

The  purpose  of  this  research  is  to  develop  a  technique  for  determining  the  three- 
dimensional  spatial  distribution  of  electron  density  in  the  ionosphere.  This  information 
is  to  be  determined  from  measurements  of  the  scattering  of  electromagnetic  waves 
entering  the  ionosphere  from  outside  the  earth's  atmosphere.  It  is  assumed  that  these 
measurements  are  to  be  ground  based,  and  that  the  source  of  the  electromagnetic 
waves  is  known.  The  availability  of  a  satellite  beacon  as  a  source  can  be  assumed.  By 
definition,  this  is  an  inverse  scattering  problem. 

As  discussed  below,  the  distinction  between  a  direct  scattering  and  an  inverse 
scattering  problem  is  based  upon  what  is  known  and  what  is  sought  in  the  problem.  In  a 
direct  scattering  problem,  the  source  of  the  fields,  the  geometry,  and  the  medium 
properties  involved  are  all  known;  the  scattered  fields  are  sought.  In  an  inverse 


scattering  problem,  the  scattered  fields  are  measured  and  thus  known.  What  is  sought 
is  a  description  of  the  sources  in  the  problem.  In  the  problem  under  consideration,  the 
ionosphere  is  probed  using  an  electromagnetic  field,  and  the  scattered  field  (or  the  sum 
of  the  incident  and  scattered  fields)  is  measured,  and  thus  known.  What  is  sought  is  the 
electron  density  as  a  function  of  three  spatial  dimensions.  Because  the  probing  wave 
interacts  with  the  density  variations  by  inducing  currents  in  the  medium  (which  may  be 
viewed  as  the  sources  of  the  scattered  field),  the  desired  information  about  the  electron 
density  profile  is  contained  in  this  induced  source  term.  Thus,  in  this  problem,  the 
fields  are  measured  and  the  (induced)  source  term  is  sought;  the  determination  of  the 
electron  density  profile  is  an  inverse  scattering  problem.  These  concepts  are  made 
clearer  in  the  mathematical  presentation  below. 

Inverse  scattering  as  a  field  has  existed  for  about  50  years.  It  originated  in  the 
quantum  mechanists'  need  to  determine  the  scattering  potential  of  an  atomic  system 
from  scattering  data.  Because  of  this  basis,  inverse  problems  have  been  thought  to 
have  inherently  ill-posed  and  nonunique  solutions  with  little  or  no  practical  application. 
This  view  has  been  due  largely  to  the  development  of  inverse  scattering  theory  from  a 
discipline  in  which  coherent  field  measurements  are  impossible,  and  in  which  it  was  the 
scattering  potential  that  was  sought  rather  than  the  source  term  in  the  wave  equation. 
Recently,  N.  N.  Bojarski  has  developed  an  Exact  Inverse  Scattering  Theory  that 
overcomes  these  limitations.  It  is  this  theory  that  forms  the  basis  for  the  research 
reported  in  this  paper.  It  is  applicable  to  quite  general  electromagnetic  and  acoustic 
problems.  As  presented  below,  its  solution  has  been  proven  unique,  well-posed,  and 
possessing  properties  which  make  it  eminantly  applicable  to  practical  applications. 

The  ionosphere  is  a  low-temperature  (but  not  "cold")  plasma  existing  in  the 
presence  of  the  earth's  magnetic  field  and  in  the  presence  of  a  complex,  temporally  and 
spatially  varying  background  of  neutral  and  ionic  components.  Literal  determination  of 
the  three-dimensional  spatial  distribution  of  electron  density  could  imply  knowledge  of 
the  position  of  each  electron  at  an  instant  of  time.  This  is  neither  possible  or  desired. 
The  purpose  for  which  a  measure  of  the  electron  density  distribution  is  wanted  is  to 
understand  and  predict  ionospheric  effects  on  electromagnetic  propagation  (e.g.,  for 
communications)  and  to  understand  ionospheric  morphology.  These  effects  depend  on 
electron  density  distributions  over  a  wide  range  of  scale  sizes.  There  is  strong  evidence 
for  this  scale  size  range  to  extend  from  less  than  1  m  to  an  appreciable  fraction  of  the 
radius  of  the  earth  (Ref  1).  For  almost  all  purposes  related  to  structural  and 
propagation  studies,  the  distances  over  which  variations  in  electron  density  are  of 
importance  are  at  least  several  meters,  and  often  many  orders  of  magnitude  larger  than 


this.  This  means  that  the  spatial  resolution  with  which  the  electron  density  needs  to  be 
determined  is  on  a  macroscopic  scale  with  respect  to  molecular  and  ionic  constituent 
interactions  with  the  probing  and  scattered  fields.  Limitation  to  such  macroscopic 
scales  is  also  justified  based  on  the  spatial  resolution  usually  available  for  field 
measurements. 

The  ability  to  view  the  electromagnetic  scattering  process- -and  thus  the  inverse 
scattering  problem--on  a  macroscopic  scale  is  of  considerable  importance.  As 
discussed  in  the  next  subsection,  the  inverse  scattering  theory  as  it  is  applied  to  this 
problem  is  based  on  the  ability  to  express  the  fundamental  physics  of  the  problem  in 
terms  of  a  linear  wave  equation  and  a  constitutive  relationship.  This  formulation  is 
quite  accurate  for  the  ionosphere  on  a  macroscopic  scale.  However,  on  a  microscopic 
scale  radiation  transport  theory  must  be  employed,  and  the  equations  are  usually 
nonlinear  and  much  more  complex.  This  is  not  the  problem  under  consideration  here. 

A  comment  on  terminology  is  in  order.  The  problem  addressed  by  this  research 
has  been  termed  the  determination  of  the  electron  density  "profile."  Historically, 
among  the  first  experimental  determinations  of  the  spatial  variation  in  electron  density 
in  the  ionosphere  were  vertical  incidence  sounder  measurements  of  the  apparent 
variation  in  density  with  height:  The  vertical  electron  density  profile.  Use  of  the  term 
"profile"  strictly  implies  a  one-dimensional  variation.  It  is  well  known  that  variations  in 
ionospheric  electron  density  in  all  three  dimensions  are  of  fundamental  importance  in 
many  problems,  and  the  research  reported  in  this  paper  is  directed  towards  obtaining 
three-dimensional  information.  However,  it  is  sometimes  easier  to  use  the  term 
"profile"  as  if  it  includes  two-  and  three-dimensional  variations  instead  of  always 
explicitly  stating  that  three-dimensional  variations  are  meant. 

1 .2  MEASURABLE  QUANTITIES,  CALCULABLE  QUANTITIES,  AND 

THE  ELECTRON  DENSITY 

On  a  macroscopic  scale,  electromagnetic  waves  interact  with  the  permeability 
and  permittivity  variations  of  the  medium.  The  measured  quantity  in  this  problem  is 
the  electromagnetic  field,  and  the  desired  quantity  is  the  electron  density  distribution. 
In  order  to  go  from  the  measured  fields  to  the  desired  electron  density  it  is  necessary  to 
have  a  set  of  equations  relating  these  quantities:  This  is  provided  by  the  inverse 
scattering  theory  presented  in  the  following  sections.  It  is  important  to  distinguish 
between  those  equations  that  represent  rigorous  electromagnetic  theory  and  those  that 
assume  a  physical  model  for  the  ionosphere. 


The  Exact  Theory  provides  equations  which  can  be  solved  for  the  total  field, 
<f>(x,w),  and  the  source  term,  p(x,w),  in  the  wave  equation  governing  propagation  in  the 
medium.  Given  that  Maxwell’s  equations  and  the  wave  equation  derived  from  them 
apply  to  the  ionosphere  (which  is  certainly  a  good  assumption  for  macroscopic 
interactions),  these  solutions  are  rigorous  and  do  not  involve  models.  The  source  term 
contains  the  effects  of  inhomogeneities  in  the  medium  in  the  form  of  induced  sources, 
due  to  the  interaction  between  the  probing  wave  and  the  medium.  In  order  to  relate  the 
sources  and  fields  to  the  macroscopic  permeability  and  permittivity  it  is  necessary  to 
have  a  constitutive  relation.  The  form  of  this  relation  chosen  in  the  derivation  below  is 
(for  the  scalar  case) 

p(x,w)  =  k2  [  n2(x,(o)-l]d>  (x,co)  (1) 

where  k  =  a  is  the  free  space  wavelength  and  n  is  the  complex  refractive  index. 

It  is  convenient  to  work  in  terms  of  the  complex  refractive  index  rather  than  the 
(equivalent)  permeability  and  permittivity.  This  form  of  the  constitutive  relation  can 
be  derived  directly  from  Maxwell's  equations  for  a  plasma  with  a  uniform  magnetic 
field  (e.g.,  Ref  2,  pp.  63-67).  No  other  assumptions  about  the  medium  are  involved  in 
arriving  at  a  constitutive  relation  of  the  form  of  Equation  1.  However,  additional 
modeling  assumptions  are  necessary  to  relate  the  refractive  index  to  the  electron 
density  and  to  the  other  physical  properties  of  the  medium.  It  will  be  shown  in 
Sect’  6  that  the  Appleton-Hartree  equation  relating  the  refractive  index  to  the 
electron  density  and  other  medium  properties  is  a  very  good  model  for  use  with  the 
methods  of  determining  the  electron  density  distribution  discussed  in  this  paper. 
Indeed,  use  of  this  model  results  in  a  particularly  simple  relationship  among  p,</>,  n,  and 
the  electron  density.  In  fact,  under  certain  conditions  that  are  usually  possible  to 
ensure  in  an  experiment,  it  can  be  shown  that  the  intermediate  steps  of  solving  for  p 
and  n  are  unnecessary.  It  is  possible  to  compute  a  quantity  directly  from  the  measured 
fields  that  is  proportional  to  the  electron  density.  This  is  discussed  in  Section  6. 

It  should  also  be  noted  that  the  inverse  scattering  theory  results  derived  in  the 
following  sections  are  applicable  whenever  propagation  in  the  ionosphere  can  be 
expressed  in  terms  of  the  wave  equation  employed  below.  Different  applications  to  the 
ionosphere  may  require  different  constitutive  relations  and/or  different  relations 
between  n  and  the  medium  properties.  However,  the  method  of  solving  for  the  source 
term  is  independent  of  the  constitutive  relationship. 
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1.3  RESULTS  OF  LITERATURE  SURVEY 

An  extensive,  computer-aided  literature  survey  was  carried  out  as  part  of  this 
research.  An  annotated  bibliography  of  references  related  to  inverse  scattering 
techniques  that  could  be  applied  to  the  problem  of  determining  the  three-dimensional 
variation  of  electron  density  in  the  ionosphere  is  included  in  Volume  II.  Based  on  this 
survey,  it  was  concluded  that  the  Exact  Inverse  Scattering  Tneory  of  N.  N.  Bojarski,  as 
applied  by  the  author  and  other  workers  in  the  field,  combined  with  the  Holographic 
Radio  Camera  technique  developed  by  the  author  (Refs  3-5),  is  best  suited  for 
addressing  the  problem  considered  in  this  research. 

The  primary  problem  with  other  inverse  scattering  techniques  is  that  they  suffer 
one  or  more  of  three  limitations:  They  are  applicable  only  to  one-dimensional 
problems;  they  are  ill-posed  in  the  presence  of  measurement  noise;  or  they  involve 
approximations  which  may  be  acceptable  for  problems  involving  some  types  of 
scatterers  (e.g.,  the  physical  optics  and  far-field  approximations  often  used  for  scatter¬ 
ing  by  conducting  bodies)  but  which  are  not  applicable  to  the  ionosphere.  Indeed,  it  can 
be  argued  with  mathematical  rigor  that  since  the  Exact  Theory  is  both  exact  and  has  a 
unique  solution  (see  Section  3),  it  is  the  only  truly  correct  approach. 

As  an  example  of  the  typical  situation  found  in  doing  the  survey,  one  other 
candidate  technique  which  has  been  applied  analytically  to  determine  the  electron 
density  profile  is  the  work  of  Jordan  and  Ahn  (see,  for  example,  Ref.  6  and  references 
therein).  Their  approach,  in  work  which  has  continued  over  the  last  decade,  is  to  use 
Kay's  and  Kay  and  Moses'  extension  of  the  Gelfand-Levitan  solution  for  one-dimensional 
inverse  problems.  In  particular,  they  are  able  to  solve  problems  in  which  the  profile  can 
be  modeled  using  an  analytic  function  with  a  small  number  of  poles  (e.g.,  one  to  three). 
Based  on  the  information  available  in  the  literature,  this  approach  suffers  from  three 
very  serious  limitations  for  application  to  the  problem  addressed  by  this  paper.  First, 
the  approach  has  not  been  demonstrated  to  be  applicable  to  problems  in  two  and  three 
dimensions.  Second,  the  behavior  of  the  approach  in  the  presence  of  measurement  noise 
has  not  been  investigated.  This  is  of  considerable  concern,  because  what  the  approach 
is  based  upon — the  Gelfand-Levitan  solution- -is  usually  considered  to  be  unstable  with 
respect  to  such  noise.  Third,  in  a  recent  paper  by  Reilly  and  Jordan  (Ref  7),  it  was 
shown  that  the  number  of  poles  required  to  approach  a  reasonable  approximation  of 
even  a  one-dimensional  ionospheric  profile  is  of  the  order  of  40,000  or  more.  The 
authors  readily  admitted  that  the  technique  was  probably  not  applicable  for  such 
problems. 
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2.  THE  EXACT  INVERSE  SCATTERING  THEORY 


2.1  THE  INTEGRAL  EQUATION  FOR  INVERSE  SCATTERING 

This  subsection  presents  a  derivation  of  N.  N.  Bojarski's  (Refs  8- 11)  "Exact 
Inverse  Scattering  Theory."  The  theory  is  derived  in  the  form  leading  to  an  integral 
equation  which  can  be  solved  numerically  in  a  stable  and  efficient  manner.  Section  4 
presents  derivations  which  parallel  that  of  this  section,  but  which  result  in  exact, 
closed-form  analytic  solutions  for  the  total  field  and  the  source  term.  This  order  of 
presentation  has  been  chosen  because  it  has  proven  to  be  the  simplest  for  newcomers  to 
the  field  of  inverse  scattering  to  understand.  It  is  important  to  note  that  although 
consideration  of  the  integral  equation  is  the  simplest  means  of  introducing  the  Exact 
Theory,  the  analytic  solutions  of  Section  4  are  by  far  the  most  significant  results 
presented  in  this  paper. 

Although  scalar  notation  is  used  in  the  derivation  of  this  subsection  and  elsewhere 
in  this  paper,  unless  otherwise  noted  all  of  the  results  presented  have  been  shown 
rigorously  to  apply  to  the  full  vector  electromagnetic  equations  for  an  anisotropic 
medium.  Specifically,  they  apply  to  the  ionosphere  in  the  presence  of  a  uniform 
magnetic  field  (see  also  the  discussions  in  Sections  1.2  and  Section  6). 

Consider  a  source  p(x)  in  a  domain  D  bounded  by  a  surface  S.  Then  the  time 
harmonic  field,  <Mx),  due  to  P(x)  is  the  solution  to  the  inhomogeneous  wave  equation 

v'vMx)  +  k ^</>(x)  =  -p(x),  x  €  D  (2) 

where  k  =  l-n/k. 

A  direct  scattering  problem  is  one  in  which  p(x)  is  known  or  specified,  and  a 
solution  for  <Mx)  is  sought.  The  inverse  scattering  problem  is  one  in  which  <f>(x)  is 
known,  and  p(x)  is  sought.  For  the  inverse  source  problem,  <f>(x)  is  measured  over  some 
surface,  and  the  object  is  to  determine  p(x).  In  general,  p(x)  =  pm(x)  +  p$(x),  where  pm 
is  due  to  interaction  with  the  medium,  and  ps  is  due  to  actual  sources.  If  n(x)  is  the 
complex  refractive  index  of  the  medium,  then 


p 


(3) 


In  most  remote  probing  problems,  Ps(x)  is  known,  and  pm(x)  *s  sought  to  yield  n(x).  This 
is  termed  the  inverse  medium  problem.  If  <Mx)  is  specified  (as  the  desired  field)  and 
p(x)  or  n(x)  is  sought  so  as  to  produce  that  <Mx),  the  problem  is  termed  an  inverse 
synthesis  problem.  This  paper  deals  with  an  inverse  medium  problem. 

Let  the  following  field  quantity,  <2>^(x),  be  defined: 


g  (x-x')vd>(x')  -<Mx')vg  (x-x')  dS' 


x-x')  ( 


(4) 


where  g(x)  is  the  free  space  Green's  function  and  the  asterisk  denotes  complex 
conjugation,  g  satisfies  Equation  2  with  p(x)  =  6  (x).  <f>H  is  in  the  form  of  the  Kirchhoff 
integral  with  g  complex  conjugated.  Note  that  if  the  Kirchhoff  integral  is  applied  to 
the  field  d»(x)  on  S  and  evaluated  at  any  point  x  inside  D,  it  is  identically  zero:  The 
Kirchhoff  integral  is  nonzero  only  for  points  outside  D.  Conversely,  is  nonzero 

only  for  points  inside  D.  Points  outside  of  D  are  associated  with  the  direct  scattering 
problem;  points  inside  D  are  of  interest  for  the  inverse  scattering  problem.  This 
topological  difference  is  the  reason  why  direct  scattering  solutions  are  mathematically 
ill-posed  when  applied  to  the  inverse  scattering  problem. 

It  should  also  be  noted  that  <f> ^  is  the  mathematical  expression  for  the  reconstruc¬ 
tion  obtained  from  a  hologram  ( <f>  in  Equation  4)  recorded  on  S.  The  relationship 
between  holography  and  inverse  scattering,  along  with  an  analysis  of  the  consequences 
for  remote  probing  and  coherent  imaging  applications,  has  been  presented  by  Stone 
(Reference  3;  see  Section  6).  6^  is,  in  general,  known  for  inverse  problems,  since  </>  is 
known  over  S.  <f>  is  measured  over  S  for  the  inverse  source  and  medium  problems,  or 
specified  over  S  for  the  inverse  synthesis  problem. 

Applying  Gauss'  theorem  to  Equation  4  converts  the  surface  integral  into  a  volume 
integral: 


dV  (g*V20-4>V2g*) 

From  Equation  2, 


(5) 
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V^>  =  -k ^</>  -P 


(6) 


and,  by  complex  conjugation  of  Equation  2  for  g, 


2  * 
V  8 


.  2  * 

=  -Kg  -5 


(7) 


Substitution  of  Equations  6  and  7  into  Equation  5  gives 


~J  dV  [s*(- 

■/ 


k2<f»-P)  -<M-k2g  -8)j 


dVW>8-g  P ) 


and,  carrying  out  the  integration  over  the  delta  function, 


</>H  =  <t>  ~  j 


dV  g *p 


Direct  scattering  theory  gives  the  result  that 


(8) 


(9) 


(f)  = 


J  dV  g  p  + 

/dVgP  +4^ 


dS(gvd>-<#>vg) 


(10) 


In  Equation  10,  the  first  integral  is  just  the  superposition  integral  over  the  sources.  The 
second  term  is  the  Kirchoff  integral,  and  is  associated  with  the  incident  field,  <£-  For 
the  inverse  scattering  problem,  can  be  assumed  to  be  known  without  loss  of 
generality  (e.g.,  it  is  the  known  probing  field  for  the  inverse  medium  case,  or  the 
specified  incident  field  in  the  inverse  synthesis  case). 

Equations  9  and  10  are  two  independent  simultaneous  equations  in  two  unknowns, 
<f>andP .  Substitution  of  Equation  10  into  Equation  9  yields 
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(11) 


=  /dvgp.y  dV  g*  p  +^. 
=  y  dV  (g-g*)  />  +<^>j 


or 


</>H(x)  =  2\  J  dV' Img(x-x')  p  (x')  +  (/>  j(x)  (12) 

where  Im  denotes  the  imaginary  part.  Equation  12  is  the  basic  equation  of  the  Exact 
Inverse  Scattering  Theory.  It  is  an  integral,  convolution  equation  for  the  single 
unknown,  p  (x).  It  can  be  solved  by  standard  deconvolution  techniques.  Quite  recently, 
Stone  has  presented  a  closed-form  solution  to  Equation  12,  based  on  a  closed-form 
solution  by  Bojarski  for  the  total  field,  <f>,  in  Equation  9  (see  Section  4). 

2.2  DETERMINATION  OF  THE  IONOSPHERIC  REFRACTIVE  INDEX 

Based  on  the  above  theory,  the  inverse  medium  problem  can  be  solved  by  the 
following  steps: 

A.  Compute  <f>^(x),  using  the  measured  field  values  in  Equation  4  (note  that  the 
surface  of  integration,  S,  is  the  measurement  surface). 

B.  Solve  Equation  12  for  p  (x),  using  $H(x)  from  A  and  the  known  incident  field, 

4>j(x). 

C.  Compute  the  total  field, 4>(x)>  from  the  direct  scattering  result,  Equation  10, 
using  p(x)  from  B. 

D.  Solve  Equation  3  for  the  desired  complex  refractive  index,  n(x),  using  p  (x) 
from  B  and  <Mx)  from  C. 

Note  that  for  the  inverse  source  problem,  only  steps  A  and  B  are  required.  However, 
for  the  inverse  medium  problem  it  is  necessary  to  carry  out  steps  C  and  D  in  addition. 
The  solution  to  the  direct  scattering  problem  (step  C)  is  a  necessary  step  in  solving  the 
inverse  medium  problem.  It  is  important  to  emphasize  that  the  computations  involved 
in  steps  A  through  C  are  all  convolution  integrals:  They  can  be  carried  out  using  fast 
Fourier  transform  techniques.  As  a  result,  computation  time  and  storage  requirements 
are  proportional  to  N  log2  N,  where  N  is  the  number  of  data  points.  Step  D  is  an 
algebraic  operation. 


The  closed-form,  analytic  solutions  of  Section  4  permit  a  closed-form,  analytic 
solution  to  the  inverse  medium  problem. 


2.3  DETERMINATION  OF  THE  IONOSPHERIC  ELECTRON  DENSITY 

When,  for  example,  the  above  four  steps  are  applied  to  fields  from  a  satellite- 
borne  beacon  which  have  propagated  through  the  ionosphere  and  are  measured  on  the 
surface  of  the  earth,  the  three-dimensional  complex  refractive  index  of  the  ionosphere 
is  obtained.  As  discussed  in  Section  1.2,  it  is  necessary  to  have  a  model  relating  the 
refractive  index  to  the  electron  density  in  order  to  obtain  the  electron  density.  As 
discussed  in  Section  6,  the  Appleton-Hartree  formula  (Reference  3,  Section  2.3.2)  is  an 
excellent  approximation  under  conditions  that  can  usually  be  readily  achieved  in 
ionospheric  measurement: 


(13) 


where,  in  the  standard  notation, 

X  =  Ne2/«omw2 
=  eB^/mai 

Yy  =  eB^./mw 

Z  =  v/<o  (14) 


with  N  being  the  (sought-after)  electron  number  density,  e  is  the  charge  and  m  is  the 
mass  of  the  electron,  is  the  permittivity  of  free  space,  co  =  2nf  is  the  angular 
frequency  of  the  probing  wave,  v  is  the  collision  frequency,  and  B  is  the  magnetic  field 
strength.  The  subscripts  L  and  T  denote  the  portion  of  the  magnetic  field  longitudinal 
and  transverse  to  the  direction  of  the  normal  to  the  propagating  wave.  In  the  absence 
of  the  magnetic  field  (Y^.=Y^=0)  and  of  collisions,  the  refractive  index  is  purely  real 
and  given  by 

n2  =  1  -  X  =  1  -  K  N/f2  (15) 
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where  K  =  e  /4ir  e  m  =  80.5  m  Hz  when  N  is  in  meters'^  and  f  is  in  Hertz.  As 

discussed  in  Section  6,  at  high  enough  frequencies  Equation  15  becomes  an  excellent 

approximation.  In  such  cases,  the  scalar  field  treatment  used  above  is  adequate,  and 

both  the  measurements  and  computations  are  simplified.  To  obtain  the  three- 

dimensional  electron  density  distribution  in  this  case,  n(x)  fr<\n  step  D  of  Section  2.2  is 

substituted  into  Equation  15,  which  is  solved  f--  N(x). 

At  lower  frequencies  the  anisotropy  of  the  ionosphere  produced  by  the  earth's 
magnetic  field  becomes  important.  Mathematically,  this  affects  two  aspects  of 
applying  the  Exact  Theory.  First,  the  above  derivation  and,  in  particular,  Equations  10 
and  12,  must  be  extended  to  include  the  full  vector  electromagnetic  fields.  The  need 
for  this  can  be  seen  immediately  from  the  form  of  Equation  13:  The  refractive  index  is 
different  for  the  components  of  the  wave  propagating  transverse  or  longitudinal  to  the 
earth's  magnetic  field  (indeed,  such  a  difference  is  one  definition  of  an  anisotropic 
medium).  Second,  there  are  two  unknown  quantities  on  the  right-hand  side  of 
Equation  13:  N,  and  v.  In  order  to  solve  for  both  (or  for  either  when  the  other  is 
unknown)  of  these  quantities,  two  linearly  independent  equations  for  n  must  be 
available.  These  are  obtained  by  measuring  the  longitudinal  and  transverse  field 
components,  carrying  out  the  solution  to  the  vector  equivalents  of  Equations  10  and  II 
and  steps  A  through  D  to  obtain  n(x)  for  the  longitudinal  and  transverse  directions,  and 
then  solving  these  two  equations  simultaneously  for  N(x)  and  v(x).  If  the  frequency  of 
the  probing  wave  is  high  enough  that  coliisional  effects  are  negligible  (but  not  so  high 
that  the  magnetic  field  effects  can  be  neglected),  then  the  two  equations  are  solved  in 
a  least-squares  sense  (in  the  presence  of  measurement  noise)  to  obtain  N(x). 

It  is  important  to  note  that  in  all  of  the  derivations  of  this  section,  and  in  all  of 
the  material  presented  in  this  paper  except  for  the  time  domain  example  discussed  in 
Section  5,  frequency  is  a  parameter.  This  means  that  measurements  and  determination 
of  n(x,w),  and  thus  of  N(x)  (which,  of  course,  is  independent  of  frequency),  can  be 
carried  out  using  a  single  frequency.  If  a  band  of  frequencies  is  available,  then  this 
additional  information  can  be  used  to  discriminate  against  noise  in  the  solution  and, 
equivalently,  to  improve  the  spatial  resolution  of  the  determination  of  N(x). 

2.4  EXACT  INVERSE  SCATTERING  THEORY  FOR  VECTOR  FIELDS 

To  derive  the  vector  analogue  of  Equations  10  and  12,  a  mixed  vector  field 
quantity  <f>H  is  defined: 
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x  (a  x  H)  -  Vg  (a*H)  +  iue  g 


’(ax  E)j 


(16) 


where  a  x  H  and  a»H  are  the  tangential  and  normal  components  of  the  magnetic  field, 
respectively,  over  the  measurement  surface  S,  and  a  x  E  is  the  tangential  component  of 
the  electric  field  over  S.  a  is  the  unit  normal  vector,  g  is  the  complex  conjugate  of 
the  free  space  Green's  function.  If  the  vector  form  of  Green's  first  theorem  is  applied 
to  Equation  16,  and  the  same  steps  as  were  carried  out  in  going  frorrr  Equation  4  to 
Equation  12  are  followed,  the  result  is  (Ref  9) 


</»H^  " 


2i  J"  dV'  ImVg(x-x')  x  3(x') 


(17) 


This  equation  should  be  compared  with  Equation  12. 

The  vector  analogues  of  the  direct  scattering  Equation  10  are  the  well-known 
Shatton-Chu  equations  for  the  E  and  H  fields  in  terms  of  3  (Reference  12,  Section  8.14). 

Use  of  Equations  16  and  17  according  to  the  four-step  solution  process  described 
in  Section  2.2  requires  a  constitutive  relation,  analogous  to  Equation  3.  The  utility  of 
the  relation  chosen  for  the  vector  case  depends,  to  some  extent,  on  the  model  that  will 
be  used  to  relate  the  medium  properties  in  the  constitutive  relation  to  the  desired 
electron  density.  As  an  illustration  that  is  often  a  very  good  model  for  the  ionosphere, 
the  magneto-ionic  theory  (which  led  to  the  Appleton- Hartree  equation  for  the  refrac¬ 
tive  index  used  above)  can  be  used.  The  vector  magneto-ionic  equations  can  be  derived 
from  either  of  two  viewpoints:  Considering  the  source  of  the  fields  to  be  conduction 
currents  (c.f.  Reference  12,  Chapter  2),  or  considering  the  source  to  be 
displacement/ polarization  currents  (c.f.  Reference  13,  Chapter  5).  Equation  17  is 
written  using  3,  which  implies  the  former  approach.  In  order  to  most  easily 
demonstrate  the  connection  to  the  refractive  index  and  electron  density  formulae  used 
above,  it  is  more  convenient  to  write  3  in  terms  of  the  displacement  current.  In 
Cartesian  coordinates,  with  z  the  direction  of  propagation  and  the  earth's  magnetic 
field  in  the  x-z  plane,  the  constitutive  relation  is  thus  (c.f.  Reference  13, 
Equations  5.31-5.32  ff) 
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Although  written  as  a  scalar  in  Equation  18,  the  refractive  index  along  the  x  and  y 

directions  does,  in  fact,  differ,  as  discussed  in  Section  2.3.  This  can  be  seen  from  the 

form  of  the  equation  for  n,  Equation  13.  The  value  of  n  depends  on  the  values  of  Yy 

and  Y,  ,  which  in  turn  depend  on  the  orientation  of  E  and  E  with  respect  to  the 
L  x  y 

magnetic  field.  To  properly  derive  the  form  of  Equation  18  thr.c  explicitly  contains  the 

transverse  (Y^  =  0  in  Equation  13)  and  longitudinal  (Yy  =  0)  values  of  n  (i.e.,  nL  and  n-j.), 

it  is  necessary  to  introduce  the  concept  of  a  polarization  matrix  for  the  fields. 

Expressions  for  Ex  and  Ey  in  terms  of  X,  Y^,  Yy,  Z,  and  the  transverse  and  longitudinal 

components  of  E  are  obtained  by  writing  E  and  E  in  terms  of  their  transverse  and 

longitudinal  components.  Rather  than  repeating  this  extremely  involved  and,  as  will  be 

shown  below,  unnecessary  derivation  here,  the  reader  is  referred  to  Reference  14, 

Sections  5.3-5.6  ff.  When  substituted  into  Equation  18,  the  result  is  two  simultaneous 

equations.  3  and  3  are  known  from  the  solution  to  Equation  17,  and  the  transverse 
x  y 

and  longitudinal  components  of  E  are  known  from  the  Shatton-Chu  equations.  The 
unknowns  are  X  and  Z,  or  the  electron  density  distribution  and  the  collision  frequency. 
These  two  simultaneous  equations  are  solved  for  these  two  unknowns. 

While  the  above  rather  involved  explanation  is  correct,  it  is  also  unnecessary  in 
practice.  It  turns  out  that  with  proper  choice  of  a  coordinate  system  the  scalar 
mathematical  treatment  of  the  previous  sections  can  be  employed.  An  orthogonal 
curvilinear  set  of  coordinates  is  defined  by  the  direction  of  the  earth's  magnetic  field 
(which  follows  a  curved  path  between  the  source  of  the  probing  field  and  the 
measurement  region),  and  longitudinal  (L)  and  transverse  (T)  coordinates  everywhere 
orthogonal  to  each  other  and  to  the  magnetic  field.  Then  longitudinal  and  transverse 
field  components  can  be  measured  on  the  surface  of  the  earth,  and  they  will  each 
satisfy  equations  of  the  form  of  Equations  10  and  12.  They  will  also  satisfy  the 
constitutive  relation  Equation  3,  with  n  given  by  n^  (Yy  =  0  in  Equation  13)  for  the 
longitudinal  field  and  n  given  by  ny  (Y^  =  0  in  Equation  13)  for  the  transverse 
component.  This  yields  directly  the  two  simultaneous  equations  for  N  and  v . 

Finally,  it  should  again  be  noted  that  proper  choice  of  the  probing  wave  frequency 
yields  a  real  refractive  index  with  N  as  the  only  unknown,  and  a  scalar  field  propagation 
problem.  This  eliminates  the  requirement  for  measuring  and  processing  vector  fields. 
This  is  discussed  in  more  detail  in  Section  6. 
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3.  UNIQUENESS 


If  the  solution  for  the  source  term  can  be  shown  to  be  unique,  then  the  solution 
for  the  electron  density  (and,  where  applicable,  the  collision  frequency)  is  unique  in 
both  the  scalar  and  vector  cases.  This  statement,  of  course,  depends  on  the  form  of  the 
constitutive  relation  and  the  model  chosen  for  the  ionosphere.  Its  truth  for  the 
magneto-ionic  model  can  easily  be  seen  from  the  discussions  in  Sections  2.3  and  2.4,  and 
the  algebraic  form  of  the  equations  involved.  The  rigorous  proof  of  the  uniqueness  of 
the  solution  for  the  source  term  is  presented  separately  below  for  the  scalar  and  vector 
equations.  In  both  cases,  it  will  be  shown  that  the  solution  is  unique  if  the  source  has 
compact  support  (i.e.,  is  identically  zero  everywhere  outside  some  finite  domain),  has 
finite  energy,  and  does  not  contain  any  nonradiating  components.  A  nonradiating  source 
is  a  source  which  produces  a  field  which  is  of  compact  support.  It  will  then  be  shown,  in 
Section  3.3,  that  nonradiating  sources  do  not  exist.  The  physical  meaning  for  the 
problem  considered  in  this  paper  of  the  compactness  requirement  on  the  support  of  the 
source  is  that  the  region  of  the  ionosphere  to  be  probed  must  be  finite. 

The  derivations  of  Section  2,  and  thus  the  uniqueness  proofs  of  this  section, 
assume  that  continuous  measurements  over  a  closed  surface  are  available  (e.g.,  that  S 
in  Equation  4  is  closed).  The  effect  of  using  discrete  (discontinuous)  measurements  over 
only  a  portion  of  the  closed  surface  is  to  affect  the  resolution  and  range  of  source 
(ionospheric)  structure  sizes  that  can  be  determined.  This  is  discussed  in  Section  6  in 
terms  of  a  band-limiting  operation  on  the  set  of  spatial  frequencies  recorded  and 
reconstructed  in  the  solution.  The  identical  spatial  band-limiting  operation  can  be 
applied  to  all  steps  of  all  of  the  uniqueness  proofs  discussed  in  this  section.  The  result 
is  that  the  solution  is  unique  for  the  measured  set  of  spatial  frequencies. 

It  should  be  noted  that  throughout  this  section  any  field  referred  to  as  being  a 
solution  of  the  wave  equation  is  also  assumed  to  satisfy  the  Sommerfeld  radiation 
condition. 
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3.1  THE  SCALAR  SOLUTION 


The  uniqueness  of  the  solution  to  Equation  12  was  first  deduced  by  Bojarski 
(Refs  9,11)  and  later  proven  more  rigorously  by  Bleistein  and  Cohen  (Ref  15).  A  simpler 
and  more  physically  understandable  proof  was  presented  by  the  author  (Ref  16).  This 
proof  is  presented  here. 

It  is  valuable  to  see  intuitively  why  the  measured  boundary  values  should  uniquely 
determine  the  source  term  in  the  wave  equation.  First  consider  a  point  scatterer 
illuminated  by  a  plane  wave  (or,  equivalently,  consider  a  point  radiator),  as  shown  in 
Figure  1.  The  scattered  field  is  a  spherical  wavefront.  Now  suppose  that  the  boundary 
values  are  to  be  measured  over  a  plane.  What  information  is  necessary  to  uniquely 
determine  the  scatterer?  Simple  geometry  shows  that  the  direction  of  the  inward 
pointing  normal  and  the  radius  of  curvature  of  the  spherical  wavefront  are  sufficient  to 
determine  the  location  of  its  center,  the  point  scatterer.  Elementary  geometry  also 
provides  that  there  is  sufficient  information  to  uniquely  determine  these  two  quantities 
if  the  height  of  the  sphere  above  the  measurement  plane  is  known  at  three  noncolinear 
points  in  the  plane.  In  electromagnetic  terms  these  heights  are  simply  the  relative 
phases  of  the  spherical  wavefront  at  the  three  points.  Thus,  measurement  of  the 
relative  phase  of  the  spherical  wavefront  at  three  noncolinear  points  in  the 
measurement  plane  uniquely  determines  the  position  of  the  point  source;  measurement 
of  the  amplitude  at  these  points  yields  the  strength  of  the  scatterer.  Now  consider  the 
physical  meaning  of  the  superposition  integral: 

0  =  j  pgdV  (19) 


This  integral  says  that  the  field  radiated  by  an  arbitrary  source  distribution  p  can  be 
represented  as  the  sum  of  the  fields,  g,  due  to  each  point  source  in  an  equivalent  point 
source  decomposition  of  P.  This,  of  course,  is  just  a  statement  of  Huygen's  principle. 
By  applying  the  reasoning  of  this  paragraph  to  each  point  source  and  the  spherical  wave 
it  generates,  it  becomes  obvious  that  there  is  enough  information  present  in  the 
scattered  field  to  uniquely  determine  the  source. 

The  above  paragraph  is  a  geometrical  interpretation  of  the  detailed  and  quite 
rigorous  mathematical  uniqueness  proof  given  by  Bleistein  and  Cohen  (Ref  15).  They 
expanded  the  source  and  the  fields  in  orthonormal  sets  of  spheroidal  wave  functions, 
and  then  showed  that  the  measured  fields  had  enough  information  to  determine  the 
coefficients  of  the  set  of  functions  describing  the  source. 
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Figure  1.  The  determination  of  a  point  scatterer,  illuminated  by  a  plane  wave, 
using  measurements  of  the  scattered  wave  at  three  noncolinear  points 
in  a  plane. 
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The  author  has  presented  (Ref  16)  a  simpler  proof  by  contradiction  (reducto  ad 
absurdum)  for  the  uniqueness  of  the  inverse  scattering  problem.  Assume  that  two 
different  source  terms,  Pj(x)  and  P2(x)  (or  two  different  refractive  index  distributions), 
produce  the  same  boundry  values,  $g(x),  on  a  measurement  surface  S.  This  assumption 
says  that  the  inverse  problem  is  nonunique.  It  follows  from  the  direct  scattering 
equation  that  the  boundary  values  due  to  the  two  sources  are 

$51  =  y*  Pj  g  dV  (20) 

$52  =  Jp2&dV  (21) 

The  difference  between  these  boundary  values  over  the  surface  is 

$SD  =  ^S1  "  ^S2 

=  y*  pt  gdv.y  P2gdv 

=  f  (PrP2) gdV  ,  (22) 

and 

$SD  =  f  P D  8  dV 

=  0  (23) 

where  Pq  is  the  difference  between  the  two  sources  ( PD  *  0  by  the  original  assumption), 
and  $SD  =  0  by  the  original  assumption.  Let  $j,  ^  and  $D  =  $r$2  denote  the  fields 
inside  S  associated  with  the  boundary  values  $S1,  $S2,  and  $SQ,  respectively.  Because 
$j  and  $2  must  both  satisfy  the  wave  equation,  Equation  2 ,  it  follows  that  $^  must  also 


satisfy  the  wave  equation.  Thus  is  a  solution  to  the  wave  equation  with  boundary 
values  which  are  0  everywhere  on  the  boundary.  By  the  same  argument  used  to  prove 
the  uniqueness  of  the  solution  to  the  wave  equation  (e.g.,  Reference  12,  Section  9.2),  it 
follows  that  must  be  identically  zero  throughout  the  bounded  volume.  But  is  the 
(identically  zero)  solution  to  the  wave  equation  with  source  term  =  ^-^2*  This 
means  Pq=  0  through  S,  which  contradicts  the  original  assumption.  Thus,  the  source  is 
uniquely  determined  by  the  boundary  values. 

3.2  THE  VECTOR  SOLUTION 

This  subsection  extends  the  author's  uniqueness  proof  for  the  scalar  case  to  the 
vector  inverse  scattering  problem. 

The  proof  for  vector  fields  and  sources  in  general,  anisotropic  media  is  identical 
to  the  scalar  proof  with  two  exceptions.  First,  the  Stratton-Chu  equations  for  the  field 
(Reference  12,  Section  8.1*0  must  be  used  in  place  of  the  scalar  superposition  integral. 
Because  the  Stratton-Chu  equations  are  linear  in  the  source  terms,  the  steps  in  the 
above  proof  still  apply.  Second,  it  is  necessary  that  the  result,  that  zero  boundary 
conditions  everywhere  over  a  closed  surface  imply  that  the  fields  are  identically  zero 
throughout  the  bounded  volume,  hold  for  fields  satisfying  the  wave  equation  for  general, 
anisotropic  media.  Fortunately,  this  again  follows  directly  from  the  proof  of  the 
uniqueness  of  the  fields  in  such  media.  Stratton  (Reference  12,  Section  9.2)  points  out 
why  this  is  so,  and  Kong  (Reference  16,  Section  7.1)  presents  an  explicit,  rigorous  proof. 

3.3  NONRADIATING  SOURCES 

In  their  proof  of  the  uniqueness  of  the  solution  to  the  inverse  scattering  problem, 
Bleistein  and  Cohen  (Ref  15)  derive  mathematically  a  class  of  sources  which  produce 
fields  that  are  identically  zero  everywhere  outside  a  finite  volume  (i.e.,  their  fields 
have  compact  support).  Since  the  fields  due  to  such  a  source  cannot  contribute  to  the 
boundary  values  which  would  be  measured  outside  the  volume,  it  would  appear  that 
these  sources  represent  a  potential  nonuniqueness  in  the  inverse  scattering  problem. 
Devaney  and  Wolf  (Ref  18)  and  Devaney  (Ref  19)  have  also  derived  and  discussed  such 
sources.  If  such  sources  did  exist  physically  (as  opposed  to  being  derivable  mathemati¬ 
cally),  they  would  pose  little  problem  for  inverse  scattering  in  the  context  of  remotely 
probing  the  ionosphere:  Nonradiating  sources,  viewed  as  induced  sources,  represent 
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nonscattering  features  of  the  medium.  However,  they  are  intellectually  bothersome. 
Fortunately,  within  the  last  year,  three  independent  proofs  that  nonradiating  sources 
cannot  exist  physically  have  been  obtained. 

The  first  proof  was  derived  by  the  author  (Ref  20)  from  the  realization  that  fields 
of  compact  support  are  physically  inconsistent  with  the  inhomr  geneous  wave  equation. 
The  proof  is  by  contradiction.  Assume  that  0  is  a  nomadiating  source  which  is  a 

proper  source  term  for  the  wave  equation.  Further,  assume  that  produces  a  field 
*^NR  *s  identically  zero  outside  some  finite  volume  V  .  By  "produces"  it  is  meant 

that  is  the  solution  to  the  inhomogeneous  wave  equation  with  source  Let  S 

be  a  surface,  enclosing  a  volume  V,  drawn  such  that  S  is  completely  outside  of  V  .  Then 
VQ  is  contained  in  V.  By  the  original  assumption,  is  zero  everywhere  on  S.  But 
then  is  a  solution  to  the  wave  equation  with  0  boundary  conditions.  It  follows,  by 
the  uniqueness  proof  for  the  solution  to  the  wave  equation,  that  <p^^=Q  throughout  V 
and,  in  particular,  throughout  V^.  But  then  (V^  +  ^^nR  =  ®  throughout  VN^,  and 
PNI^=0  throughout  V^.  This  contradicts  the  original  assumption.  This  proof  would 
fail  if  the  particular  frequency  involved  happened  to  be  an  eigenfrequency  for  the 
volume  V,  except  that  the  conditions  must  be  true  simultaneously  for  all  possible 
volumes  containing  V^,  and  ail  surfaces  surrounding  V. 

The  second  proof  was  obtained  by  the  author  during  the  course  of  the  research 

reported  in  this  paper.  The  basic  form  of  the  exact  analytic  solution  for  the  fields  in 

the  inverse  scattering  problem  obtained  by  Bojarski  (see  Section  4)  is  an  expression 

which  shows  that  the  spatial  Fourier  transform  of  the  field,  (k,uj),  is  given  by  e(k,u,v), 
2  2  2 

evaluated  at  v  =  u  /k  (on  the  Ewald  sphere,  in  the  terminology  of  Section  4).  It  is 
also  shown  in  Section  4  th*.t  the  Fourier  transform  of  the  field  is  identically  zero  for 
v  *  u  /k  .  This  means  that  the  spatial  Fourier  transform  of  the  field  ir  of  compact 
support.  It  is  a  well-known  property  of  the  Fourier  transform  (c.f.  Reference  21, 
pp.  143  ff)  that  if  the  transform  of  a  function  has  compact  suppoit,  then  the  function 
cannot  have  compact  support.  It  follows  that  fields  which  are  proper  solutions  of  the 
inhomogeneous  wave  equation  cannot  be  of  compact  support.  Since  the  compactness  of 
the  support  of  the  fields  produced  by  a  nonradiating  source  is  a  necessary  condition  for 
the  existence  of  such  sources,  this  means  nonradiating  sources  cannot  exist. 

The  third  proof  was  suggested  by  N.  N.  Bojarski  (Ref  22).  In  a  recent  paper, 
Porter  and  Devaney  (Ref  23)  prove  that  the  solution  to  the  inverse  source  problem 
which  possesses  minimum  source  energy  is  the  solution  in  which  any  nonradiating  part  is 
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zero.  Bojarski  has  pointed  out  that  the  physical  meaning  of  Lenz's  law  is  that  currents 
and  charges  always  assume  the  minimum  energy  distribution.  Hence,  any  nonradiating 
source  contribution  is  zero. 

In  summary,  the  solution  to  the  inverse  scattering  problem  is  unique,  in  both  the 
scalar  and  the  vector  cases. 


22 


4.  EXACT,  CLOSED-FORM,  ANALYTIC  SOLUTIONS  FOR  THE 
TOTAL  FIELD  AND  THE  SOURCE 

4.1  THE  SOLUTION  FOR  THE  FIELD  (Refs  24,25) 

In  the  derivation  of  Section  2.1,  Equations  9  and  10  were  arrived  at  as  two 
simultaneous  equations  in  two  unknowns:  The  sources,  p,  and  the  total  field,  0.  It  was 
chosen  to  eliminate  p  and  solve  for  p.  Bojarski  has  shown  (Ref  25)  that  by  starting  with 
a  slightly  different  form  of  the  initial  equations  and  eliminating  p  instead  of  p,  it  is 
possible  to  obtain  a  closed-form,  analytic  solution  for  the  total  field,  <t> . 

The  derivation  begins  with  a  somewhat  more  general  form  of  the  wave  equation 

V2p(x,w)  +  jw2/c2(x,oj)j0(x,w)  =  -p(x,u>)  (24) 

where  w  is  the  temporal  frequency  (which  will  often  be  suppressed  in  this  section),  and  c 
is  a  velocity  which  varies  as  a  function  of  both  x  and  w.  If  the  constitutive  equation  is 
written  as 


Pt(x,«,v)  =  ||w2/c2(x,w)j  -  «^/v2^0+  p(x,w)  (25) 

where  v  is  an  arbitrarily  chosen  constant  "free  space"  reference  velocity,  then 
Equation  24  can  be  written  as 

V2p(x,w)  +  |w2/v2j0(x,u>)  =  -  Pt(x,w,v)  (26) 

Pt  is  introduced  to  permit  writing  Equation  30  with  a  constant  coefficient  in  front  of  0. 
If  G  is  the  Green's  function  associated  with  the  constant  reference  velocity,  v,  then  G 
satisfies  Equation  26  with  a  delta  function  source.  Note  that  G  satisfies  Equation  26, 
and  not  Equation  24. 

Let  a  quantity  6  be  defined  as 
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(27) 


0  = 


f 


dS  .  (GrV<p-0VGr) 


where  Gr  =  ReG,  which  also  satisfies  Equation  26  with  a  delta  function  source.  6  is 
analogous  to  0j_j  of  Equation  4,  and  has  many  of  the  same  properties.  In  particular,  9 
depends  only  on  the  measured  data,  and  is  thus  known.  If  the  same  steps  that  were 
carried  out  in  going  from  Equation  4  to  Equation  9  are  applied  to  Equation  27,  the 
following  integral  equation  is  obtained: 


e 


dVGrPt 


(28) 


Equations  9  and  28  are  quite  similar. 

Let  a  new  source  term  py  be  introduced  which  has  support  only  over  the  volume 
of  integration  V: 


Py  (x,u>,v)  =  Pt  (x,0J,v)  ,  x  e  V 

0  ,  x  €  V 


(29) 


Then  Equation  28  can  be  written  as 


r 

0(x, u>,v)  =  0(x,w)  -  /  Gr(x|x',to,v)  p^(x',w,v)dnx' 

«/— OO 


(30) 


In  Cartesian  coordinates  Gr  is  a  difference  kernal  and  the  integral  in  Equation  30 
becomes  a  convolution.  Using  superscript  tildes  to  denote  the  n-dimensional  spatial 
Fourier  transform  of  quantities,  the  spatial  Fourier  transform  of  Equation  30  becomes 

<T(k,w,v)  =  0(k,w)  -  Gr(k,a;,v)p^(k,w,v)  (31) 


where  k  is  the  transform  variable.  Note  that  0  and  its  spatial  Fourier  transform  cannot 
depend  on  the  arbitrary  reference  velocity  v. 

The  spatial  Fourier  transforms  of  the  real  and  imaginary  part  of  the  Green's 
function  are 
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(32) 


Gr(k,^,v)=  P  l/(k2-«?/v2) 
l/(k2-o»2/v2) 


kW/v2 

k2=  w2/v2 


G^k.w.v)  =  U/2k)|&(k-ty)  -  5  (  k+$  )J 


where  k  =  Ikl  and  P  denotes  the  principal  value.  Note  that  the  support  of  G  is 

2  2  2  ~~  r 
everywhere  except  on  the  sphere  k  =  co  /v  ,  and  the  support  of  G.  is  only  on  the  sphere 

k2  =  Jh2.  This  sphere  is  termed  the  Ewald  sphere.  Note  also  that  the  spatial  Fourier 

transform  of  G  is  invariant  with  respect  to  the  dimensionality  of  the  space.  Because 

the  support  of  Gr  is  everywhere  except  on  the  Ewald  sphere, 


Gr(k^u,v) 


v=Wk  =  Gr(k,ai,a,/k)  =  0 


(33) 


2  2  2 

If  Equation  31  is  evaluated  for  v  =  w  / k  (on  the  Ewald  sphere),  then,  by 
Equation  33, 


(jf  =  <£(k,<u) 


(34) 


Taking  the  inverse  spatial  Fourier  transform  of  Equation  34  yields 

OO 

=  (l/2rr )nJ  e'‘-*-ff(k,to,w/k)  dnk  (35) 

—  OO 

This  is  an  exact,  closed-form  analytic  solution  for  the  total  field,  <f> .  Note  that  <!>  does 
not  depend  on  the  arbitrary  reference  velocity,  v,  used  in  the  Green's  function,  since 
must  satisfy  Equation  24,  which  is  independent  of  v. 

Previous  solutions  to  the  inverse  scattering  problem  have  yielded  expressions  for 
the  source  which  were  dependent  on  the  velocity  used  in  the  Green's  function.  Without 
additional  information,  determination  of  the  source  off  the  Ewald  sphere  surface  is 
precluded.  The  Jthor  (Ref  26)  has  shown  that  the  requirement  that  the  source  be  of 
compact  (finite)  support  is  both  necessary  and  sufficient  additional  information. 
However,  the  solution  given  by  Equation  35  achieves  the  same  effect  in  an  elegant 
fashion.  As  the  free  parameter  v  is  varied,  the  Ewald  sphere  on  which  ff  is  evaluated 
sweeps  throughout  spatial  Fourier  transform  space  to  obtain  the  complete  spatial 
reconstruction  of  </>(x). 
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A  more  complete  derivation  of  these  results  is  given  in  Bojarski  (Ref  25). 


*.2  THE  SOLUTION  FOR  THE  SOURCE 

This  result  is  perhaps  the  most  important  result  obtained  by  the  author  in  the 
course  of  the  research  reported  in  this  paper. 

A  major  advantage  of  solving  the  inverse  scattering  problem  using  numerical 
solution  of  the  integral  equation,  Equation  12,  as  discussed  in  Section  2,  is  that  fast 
Fourier  transform  (FFT)  techniques  can  be  employed.  This  makes  treatment  of  three- 
dimensional,  real-world  problems  well  within  standard  computer  capabilities.  However, 
numerical  evaluation  of  Equation  35  is  not  necessarily  as  efficient,  because  of  the 
problems  associated  with  evaluating  o'.  Furthermore,  although  Bojarski  (Ref  25)  has 
shown  that  by  evaluating  the  Laplacean  of  Equation  35  in  the  Fourier  domain  (using  the 
Fourier  differentiation  theorem)  an  expression  for  the  source  can  be  obtained,  it  would 
be  preferable  to  have  an  expression  for  the  source  directly  in  terms  of  the  measured 
data. 

The  author  has  derived  two  results  presented  here.  First,  an  expression  for  (T  in 
terms  of  Fourier  transform  operations  on  the  measured  fields  is  presented.  This  makes 
numerical  evaluation  of  the  solution  for  the  fields  using  FFT's  possible.  Second,  an 
expression  for  the  source  in  terms  only  of  (T (and  thus  the  measured  data)  is  obtained. 

Bojarski  (Ref  25)  has  pointed  out  that  since  on  the  Ewald  sphere  shell  Gr  =  0,  for 
this  case 


0 


dSVG  b 
-  r 


(36) 


This  avoids  the  requirement  to  evaluate  Vp  over  the  measurement  surface,  which  is 
quite  desirable  when  the  measured  field  contains  noise.  Consider  the  case  where  the 
measurement  surface  is  the  x-y  plane  (more  general,  nonplanar  surfaces  can  be 
accommodated  using  the  same  formulation  by  adjusting  the  phase  at  each  measurement 
point  such  that  the  measured  data  is  the  same  as  if  it  had  been  measured  over  a  plane). 
Then  Equation  36  becomes  (cj  and  v  are  temporarily  suppressed) 


0(x,y,z) 


dx'dy'  0(x',y',O)  d/dz  Gr(x-x',y-y',z) 


(37) 
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Evaluation  of  Equation  37  using  numerical  integration  requires  spatially  sampling  <f>  on 
the  measurement  surface  at  an  interval  small  compared  to  oscillations  in  Gr,  and  thus 
small  compared  to  a  wavelength.  This  is  impossible  in  most  real-world  situations; 
numerical  integration  is  also  slow.  These  limitations  can  be  avoided  by  evaluating  0  in 
the  Fourier  transform  domain,  using  FFTs.  This  requires  an  analytical  expression  for 
the  two-dimensional  spatial  Fourier  transform  of  the  norma!  derivative  of  Gr*  Noting 
that  2Gf  =  G+G*,  this  reduces  to  evaluating  1/2  3/3 z  (G+G*).  Harrington  (Ref  27) 
derives  the  following  expression  for  G*: 

d(27ri>x)d(27rj/y) 

(38) 

2  o  2Yi 

where  q=  co/v  and  r  =  (x  +  y  +  z  )  .  Taking  the  derivatives  of  both  sides  of 
Equation  38  with  respect  to  z  and  applying  the  definition  of  the  Fourier  transform 
yields 

-vja/azCj  =^,{3/8z(w)| 

=  <-l/27T)e-iz 
Sherman  (Ref  32)gives 


f  denotes  the  two-dimensional  spatial  Fourier  transform  from  x,y  to  v  ,  v  • 
xy  x  y 

Equation  37  is  a  convolution  in  x  and  y.  It  follows,  using  Equations  39  and  40  and 
the  convolution  theorem  that 


27 


r 


where 


0(x,y,z,a>,v)  =  /‘y  |  /xy|<^(x,y,0,a>)^  P(i>x,  tyz,  w  ,v) 


(41) 
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This  expression  is  strictly  valid  for  v  =  w/k.  A  similar  expression,  involving  an 
additional  term  of  the  same  form  as  in  Equation  41,  has  been  derived  for  v^cu/k. 

The  importance  of  Equation  41  is  that  it  permits  evaluating  0  (and  thus<H  using 
FFTs.  Indeed,  the  algorithm  is  quite  simple: 

1.  Compute  the  two-dimensional  spatial  Fourier  transform  of  the  measured  field 
data,  0(x,y,O),  using  an  FFT. 

2.  Multiply,  in  the  spatial  frequency  domain,  by  the  "propagator"  function,  P, 
given  in  Equation  42. 

3.  Compute  the  two-dimensional  inverse  spatial  Fourier  transform  of  the 
product,  using  an  FFT. 

This  method  of  evaluating  0  is  analogous  to  the  Fourier  domain  method  often  used  for 
reconstructing  holograms  (c.f.  Section  6),  and  for  evaluating  the  Kirchhoff  integral 
(Ref  28), 

Given  an  efficient  means  of  computing  0 ,  the  evaluation  of  (f>  becomes  equally 
efficient.  To  obtain  an  equally  efficient  expression  for  the  source  term,  py,  substitute 
Equation  32  into  Equation  31.  Requiring  that  k^tw/v,  this  yields 


0(k,w,v)  -  $(k,a>)  -  (k2-o^/v2)"^  £T(k,ci),v) 


(43) 


Substituting  Equation  34  into  Equation  43  and  carrying  out  some  algebra  gives 


p^(k,cti,v)  r  [<r(k,a>,w/k)  -  <T(k,a>,v)]  (k2-w2/v2) 


(44) 
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Equation  44  is  an  exact,  closed-form  analytic  solution  for  the  source  term  in  the  inverse 
scattering  problem.  Physically,  it  says  that  p y  can  be  written  as  the  difference 
between  the  8  field  on  the  Ewald  sphere  and  the  8  field  off  the  Ewald  sphere.  Although 
they  will  not  be  repeated  here,  all  of  the  results  in  the  other  sections  of  this  paper 
concerning  the  nonexistence  of  nonradiating  sources,  the  effects  of  incomplete 
measurements  and  noise,  the  well-posedness  of  the  solution,  .tc.,  can  be  deduced  from 
this  analytic  solution. 


5.  SIMULATION 


At  the  start  of  this  research,  one  goal  was  to  carry  out  a  reasonably  realistic  two- 
or  three-dimensional  numerical  or  analytical  simulation  of  the  determination  of  the 
refractive  index  profile  using  the  inverse  scattering  techniques  of  Section  2.  This  goal 
was  arrived  at  based  on  the  assumption  that  analytical  and/or  numerical  methods  for 
computing  the  necessary  directly  scattered  fields  were  available.  Such  methods  are  a 
necessary  prerequisite  for  any  simulation,  since  the  fields  produced  by  the  interaction 
of  a  probing  wave  as  it  propagates  through  the  ionosphere  represent  the  (simulated) 
"measured"  data  used  as  input  to  the  inverse  scattering  solution.  As  it  turned  out,  the 
research  tested  this  hypothesis,  and  found  that  such  methods  did  not  really  exist.  This, 
coupled  with  the  resources  devoted  to  the  very  important  and  unforeseen  analytical 
results  presented  in  Section  4,  resulted  in  only  partial  accomplishment  of  the  simulation 
goal.  However,  four  substantial  and  very  useful  results  were  obtained: 

1.  The  numerical  solution  to  the  inverse  scattering  integral  equation  of  Sec¬ 
tion  2  was  implemented  for  the  two-dimensional  case.  In  the  process,  several 
preliminary  results  concerning  a  very  broad  class  of  deconvolution  problems 
were  obtained. 

2.  Both  numerical  and  analytical  simulations  employing  the  Epstein  profile  were 
attempted.  It  was  finally  realized  that  this  approach  is  not  applicable  to  a 
trans-ionospheric  experiment,  and  the  reasons  for  this  were  determined. 
Based  on  rather  extensive  discussions  with  workers  in  the  ionospheric 
propagation  field,  the  author  believes  that  this  point  is  not  well  known. 

3.  A  method  for  calculating  the  required  directly  scattered  data,  to  use  as  input 
to  the  inverse  scattering  simulation,  was  developed.  Its  implementation  was 
completed,  but  its  testing  was  not. 

4.  To  demonstrate  the  features  that  should  be  obtainable  with  the  inverse 
scattering  technique  of  Section  2,  a  one-dimensional  simulation  was  imple¬ 
mented.  This  employed  a  time  domain  direct  and  inverse  scattering  theory 
that  can  be  viewed  as  analogous  to  the  theory  presented  in  Section  2.  Results 


were  obtained  for  quite  realistic  one-dimensional  ionospheric  proxies.  The 
ability  to  obtain  accurate  determination  of  the  profile  between  locet  maxima 
and  above  the  absolute  maximum  of  the  electron  density,  using  only  data  that 
could  be  obtained  with  a  single,  vertical  incidence  sounder  operating  on  the 
ground,  was  demonstrated.  This  has  the  potential  'or  direct  application  to 
the  processing  of  existing  ionospheric  sounder  data 

The  above  results  are  reported  in  Sections  4.1  through  4.4,  respectively.  Finally, 
it  is  pointed  out  in  Section  4.5  that  a  recent  result  of  the  author's  (presented  in 
Section  6)  shows  that,  under  conditions  that  would  probably  be  employed  in  a  practical 
ionospheric  experiment,  results  from  previous  numerical  and  experimental  work 
achieved  the  currently  desired  simulation  goal  in  a  manner  not  previously  recognized. 

5.1  NUMERICAL  SOLUTION  TO  THE  INVERSE  SCATTERING  PROBLEM 

It  was  shown  in  Sections  2.2  and  2.3  that  solving  the  inverse  scattering  problem  to 
obtain  the  electron  density  distribution  consists  of  five  steps.  The  first  step  is  to  solve 
the  basic  integral  equation,  Equation  12,  for  the  source,  P(x).  The  subsequent  four  steps 
are  straightforward  to  implement;  this  first  step  is  not. 

In  two  dimensions,  Equation  12  is  written  in  discrete  form  as  follows: 

N-l  M-l  5) 

m)  =  ^2  h(^-p,m-q)p(p,q) 

p=-N+l  q=-M+l 

where  h  is  Img,  and  discrete  indices  have  been  substituted  for  Cartesian  variables, 
implying  (usually  uniform)  sampled,  discrete  values.  Thus,  <f>H(C,m)  is  the  complex 
value  of  <#>^j  at  x  g,ym  within  the  N  by  M  set  of  x,y  points  at  which  <f>^  is  known.  and 
p  in  Equation  45  are  two-dimensional  N  by  M  matrices,  h  is  a  block,  or  vector,  Toeplitz 
matrix.  A  comprehensive  discussion  of  Toeplitz  matrices  is  given  in  Reference  29.  A 
Toeplitz  matrix  is  a  matrix  in  which  elements  along  each  left-to-right  diagonal  are  the 
same.  A  block  Toeplitz  matrix  consists  of  submatrix  blocks  which  are  Toeplitz, 
arranged  in  turn  in  Toeplitz  fashion  (as  shown  in  Equations  46  through  49,  below).  It  can 
be  shown  that  a  matrix  will  have  this  property  if  matrix  elements  for  which  the 
difference  in  the  row,  column  indices  is  the  same  have  the  same  value.  Because  the 
Green's  function  depends  only  on  the  difference  of  the  coordinates,  the  discrete  matrix 
representation  of  Img  has  this  property.  Indeed,  any  two-  or  three-dimensional 
convolution  relationship  leads  to  an  equation  of  the  form  of  Equation  45.  Because  of 
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this,  several  of  the  results  obtained  below  can  be  applied  to  any  linear  system  in  which 
the  output  is  expressable  in  terms  of  the  convolution  of  the  input  with  an  "impulse 
response"  function. 

Equation  45  can  be  written  in  matrix  notation  as  follows  (a  single  underline 
denotes  a  column  vector;  a  double  underline  denotes  a  matrix). 


±H  = 


(46) 


where 


h  =  (hP*^)  0<p,q<N-l,N>l 


hr  =  (hr  )  0<p,q<M-l,M>l 

—  p-q 


(47) 

(48) 


and 

h^  =  h(r,s)  (49) 

In  forming  the  MN  dimensional  column  vectors  ^  and  £,  the  first  N  elements  of 
for  instance,  are  the  elements  of  the  first  row  of  (f>H(C,m),  the  next  N  elements  are 
from  the  second  row,  etc.,  for  M  rows. 

Equation  46  can  be  solved  for  £  by  inverting  the  matrix  h.  This  procedure  was 
implemented,  using  a  standard,  very  well  written  double  precision  complex  matrix 
inversion  routine  which  employed  Gaussian  elimination  with  pivoting  and  equilibration. 
The  resulting  routine  was  tested  using  various  combinations  of  delta  function  source 
terms,  and  the  associated  values  of  4>^  computed  by  direct  evaluation  of  Equation  45. 
The  solutions  for  p  were  accurate  to  full  single  precision  machine  accuracy  when  single 
precision  values  of  0^  were  employed  and  no  simulated  measurement  noise  was  added 
to</>^.  The  fact  that  a  double  precision  matrix  inversion  routine  is  required  to  yield 
single  precision  accuracy  is  due  to  two  causes.  First,  Gaussian  elimination  with 
pivoting,  even  with  equilibration,  suffers  from  a  significant  loss  of  precision  unless  the 
matrix  being  inverted  is  very  well  conditioned.  Although  this  has  been  known  in 
numerical  analysis  circles  for  many  years  (c.f.  Ref  30),  very  few  computer  routines 
consider  the  conditioning  of  the  matrix,  and  almost  none  apply  the  required  iterative 


improvement  techniques  (Ref  30).  This  is  true  in  spite  of  the  fact  that  a  properly 
written  routine  which  employs  iterative  improvement  can  readily  yield  single  precision 
accuracy  for  the  inverse  using  only  single  precision  elimination,  pivoting,  and  other 
computations.  Only  after  considerble  effort,  and  only  as  this  paper  is  being  written, 
was  the  author  able  to  obtain  a  proper  inversion  routi  ie  employing  iterative 
improvement. 

The  inherent  ill-posedness  of  Equation  45  when  h  is  derived  from  the  free  space 
Green's  function  (or  its  imaginary  part)  is  well  known.  It  is  a  result  of  the  many  powers 
of  10  variation  in  the  numerical  values  of  h  which  occur  when  propagation  distances 
larger  than  a  few  wavelengths  are  involved.  Fortunately,  the  solution  to  Equation  45 
has  been  shown  by  the  author,  both  numerically  and  analytically,  to  be  rendered  well- 
posed  when  certain  a  priori  information  is  employed  in  solving  for  it.  This  is  true  even 
in  the  presence  of  measurement  noise  in  <£j_|  (see  Section  6).  The  primary  piece  of 
information  necessary  to  obtain  this  well-posedness  is  the  fact  that  the  source  is  of 
compact  support.  Since  solutions  over  finite  volumes  utilizing  data  measured  over 
finite  apertures  are  all  that  are  sought  in  any  real-world  situation,  this  information  is 
always  available. 

To  summarize,  the  ability  to  numerically  solve  Equation  12  for  the  source  term 
was  demonstrated.  The  behavior  in  the  presence  of  noise  was  consistent  with  previous 
results,  as  discussed  in  Section  6.  Extensive  experimentation  with  the  behavior  of  the 
discrete  form  of  Equation  12  was  limited  by  the  unavailability  of  a  proper  matrix 
inversion  routine.  Such  a  routine  has  now  been  obtained,  and  is  available  for  future 
work. 

In  the  process  of  examining  the  numerical  behavior  of  equations  of  the  form  of 
Equation  45,  several  additional  but  preliminary  results  were  obtained.  Because  Equa¬ 
tion  45  has  broad  application  to  almost  any  linear  physical  system,  these  results  deserve 
pursuit  in  future  work  (although  not  necessarily  as  a  part  of  additional  research  on  the 
ionospheric  profiling  problem). 

As  a  part  of  the  process  of  understanding  the  numerical  properties  of  Equations  12 
and  45,  the  author  carried  out  in-depth  research  on  the  inversion  of  block  Toeplitz 
systems  in  the  presence  of  noise.  The  following  paragraphs  summarize  this  work. 

Many  papers  on  the  inversion  of  Toeplitz  systems  have  appeared.  Howeve.  >nly  a 
relatively  few  of  these  deal  with  block  Toeplitz  systems,  and  are  thus  applicable  to 
two-  and  three-dimensional  problems.  The  inversion  of  block  Toeplitz  matrices  is 
accomplished  via  one  of  two  general  approaches:  Approximation  via  a  circulant  matrix 
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and  inversion  using  fast  Fourier  transform  (FFT)  methods,  or  direct  inversion  taking 
advantage  of  the  symmetry.  The  approach  chosen  has  very  important  implications  for 
both  efficiency  and  accuracy. 

Real-world  measurements  can  involve  anywhere  from  32  by  32  data  arrays  to 

arrays  of  several  hundred  square.  The  inversion  time  depends  on  the  product  of  the 

dimensions  (or  the  total  number  of  elements)  of  the  Toeplitz  matrix.  Let  this  quantity 

be  L  (=1024  for  a  32  x  32  matrix).  Standard  matrix  inversion  techniques  require  times 

proportional  to  with,  for  example,  one  minute  of  execution  required  for  L  =  256,  but 

over  an  hour  required  for  a  32  x  32  matrix  on  a  VAX  11/780.  (These  times  are  based  on 

the  requirement  to  use  double  precision,  as  discussed  above.  A  factor  of  ten  or  more 

savings  should  be  realized  by  going  to  single  precision  methods  using  iterative 

improvement.)  Algorithms  which  take  advantage  of  the  symmetry  of  a  Toeplitz  matrix 

2 

have  reduced  this  over  the  last  decade  to  a  constant  times  L  .  The  reader  can  assess 
this  progression  and  the  current  status  of  such  algorithms  by  consulting  References  31 
through  42.  Many  of  the  published  algorithms  have  not  been  implemented,  or  their 
implementations  are  not  applicable  to  the  solution  of  Equation  45.  In  short,  no  "direct" 
inversion  algorithm  which  exploits  the  structure  of  a  Toeplitz  matrix  and  which  is 
applicable  to  block  Toeplitz  forms  could  be  found.  Furthermore,  the  behavior  of  these 
algorithms  in  the  presence  of  ill-conditioning  and  noise  is  largely  unknown. 

One  very  promising  development  has  occurred  in  this  field  during  the  period  of 

this  research.  Three  groups,  working  independently,  have  developed  algorithms  which 
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claim  to  invert  Toeplitz  matrices  in  times  of  order  L  log  L  (Refs  43-46).  Unfortun¬ 
ately,  discussions  with  at  least  one  author  from  each  of  the  three  groups  have 
established  that  as  of  duly  1981  none  of  the  algorithms  has  been  implemented  for 
Toeplitz  systems,  let  alone  for  block  Toeplitz  matrices.  It  is  anticipated  that  a  block 
Toeplitz  routine  may  be  available  circa  January  1982. 

As  this  paper  was  being  written,  the  author  was  able  to  obtain  a  routine  written 
by  S.  Wood  (Ref  47)  (who  did  her  Ph.D.  under  M.  Morf--c.f.  Ref  43)  which  provides  a 
partial  implementation  of  these  methods  for  block  Toeplitz  matrices.  There  has  been 
no  opportunity  to  examine  this  routine  in  detail. 

A  circulant  matrix  is  a  matrix  in  which  the  elements  of  each  row  are  the  elements 
of  the  previous  row,  right-circularly-shifted  one  position.  Circulant  matrices  can  be 
inverted  in  order  L  log2  L  operations,  because  their  eigenvalues  can  be  obtained  by 
computing  the  discrete  Fourier  transform  of  a  single  row  or  column.  Toeplitz  matrices 


play  a  fundamental  role  in  digital  image  processing  and  restoration.  However,  the 
minimum  size  matrix  of  interest  for  such  applications  tends  to  be  of  order  256  x  256,  or 
larger.  This  has  led  to  extensive  development  of  fast  approximate  techniques  for 
inversion,  based  on  the  approximation  of  the  Toeplitz  matrix  by  a  related  circulant 
matrix.  Reference  43  (Chapters  7  and  8)  and  Reference  44  (Chapter  5)  review  such 
techniques  in  detail,  as  well  as  providing  over  two  dozen  references  used  by  the  author. 
The  fundamental  problem  in  applying  these  techniques  to  the  inverse  scattering  problem 
is  that,  although  the  Fourier  transform  operations  tend  to  improve  the  inversion  over 
other  methods  in  the  presence  of  ill-conditioning,  the  approximation  of  the  original 
Toeplitz  matrix  introduces  significant  errors.  This  is  in  part  due  to  the  rather  heurestic 
methods  used  to  arrive  at  the  circulant  approximation  of  the  Toeplitz  matrix  (an 
infinity  of  such  approximations  exist  for  any  particular  Toeplitz  matrix).  There  are  two 
interesting  exceptions  to  this,  however. 

Gray  showed  (Ref  50)  that  the  eigenvalues  of  a  finite  Toeplitz  matrix  and  certain 
of  its  circulant  approximations  are  asymptotically  equal.  In  a  paper  that  has  apparently 
been  almost  totally  overlooked  by  the  community  dealing  with  Toeplitz  matrices, 
Ekstrom  (Ref  51)  combined  this  choice  of  circulant  approximation  with  an  iterative 
improvement  method  (such  as  discussed  above)  and  applied  the  result  to  the  inversion  of 
block  Toeplitz  matrices.  The  result  has  an  execution  speed  of  the  order  of  the  FFT 
within  each  iteration,  and  appears  to  be  capable  of  handling  ill-conditioned  matrices. 
There  was  insufficient  time  remaining  to  implement  Ekstrom's  algorithm  by  the  time 
the  basis  for  its  operation  and  its  significance  were  understood  by  the  author.  However, 
this  approach  should  receive  high  priority  in  future  work. 

Equation  45  is  a  two-dimensional  convolution  equation.  It  can  be  solved  using  FFT 
techniques.  Let  f  and.f  denote  the  direct  and  inverse  discrete  Fourier  transform  of 
a  quantity,  respectively.  Then  the  solution  to  Equation  45  can  be  computed  as 

p  = /’1^H]/7[h]|  (50) 

The  problems  of  using  Equation  50  when  (j> ^  contains  noise  and  when  the  Fourier 
operations  are  indeed  discrete  (and  therefore  approximate)  are  discussed  in  detail  in 
Reference  48  (Sections  6.1  and  6.2)  and  Reference  49  (Section  5.3.1  ff).  It  has  been 
shown  by  this  author  that  the  approach  of  Equation  50  can  be  made  well-posed  by 
applying  constraints  using  a  priori  information  that  is  always  available  (Refs  52,53). 
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However,  a  general  algorithm  for  applying  these  constraints  is  not  yet  available.  Hunt 
(Ref  54;  Ref  48,  pp.  148-50  ff;  Ref  49,  Section  5.6)  has  developed  an  algorithm  for 
image  restoration  which  carries  out  the  processing  of  Equation  50  using  a  constrained 
least-squares  method  in  which  the  constraint  utilizes  information  on  the  statistics  of 
the  noise  in  the  data  which  would  be  available  in  any  practical  inverse  scattering 
experiment.  In  one  sense,  it  can  be  proven  that  this  method  yields  the  best  solution 
obtainable  in  the  presence  of  noise.  There  is  also  a  close  connection  between  the 
method  of  applying  the  constraint  in  Hunt's  algorithm  and  the  method  employed  by  the 
author  in  proving  the  above-mentioned  results  on  well-posedness.  This  approach  should 
also  be  seriously  examined  in  future  work. 

5.2  INVESTIGATIONS  OF  DIRECT  SCATTERING  SIMULATION  USING  EPSTEIN 

AND  RELATED  PROFILES 

A  significant  amount  of  work  has  been  carried  out  modeling  ionospheric  propaga¬ 
tion  using  Epstein  and  related  profiles.  These  are  profiles  which  vary  in  one  dimension 
according  to  an  analytical  formula.  The  formula  for  the  profile  appears  in  the 
coefficient  multiplying  the  field  in  the  wave  equation.  By  proper  choice  of  the  analytic 
form  of  the  profile,  the  wave  equation  takes  the  form  of  an  equation  that  can  be  solved 
analytically  in  terms  of  known  functions.  Usually,  a  profile  involving  the  hyperbolic 
trigonometric  functions  is  used,  and  the  wave  equation  is  transformed  into  the 
hypergeometric  equation.  Examples  of  this  approach  which  were  examined  in  depth 
during  this  research  are  contained  in  References  55  through  69.  Books  by  Brekhovskikh 
(Ref  70)  and  Wait  (Ref  71)  treat  this  topic  in  detail.  This  approach  appears,  at  first 
glance,  to  be  ideal  for  computing  the  direct  scattering  data  to  be  used  as  input  to  an 
inverse  scattering  simulation.  It  gives  an  analytic  solution.  There  are  two  problems, 
however:  The  approach  is  one  dimensional,  and  it  applies  to  reflection  problems  only. 

A  two-dimensional  problem  of  interest  can  be  treated  using  the  same  mathe¬ 
matics  by  employing  a  point  source  as  the  illumination.  Although  most  of  the  one¬ 
dimensional  problems  treated  in  the  literature  assume  plane  wave  illumination  incident 
normally  on  the  medium,  solutions  for  a  plane  wave  at  an  arbitrary  angle  of  incidence 
are  available.  The  point  source  can  then  be  expanded  into  an  angular  spectrum  of  plane 
waves,  and  the  arbitrary  angle  of  incidence  solution  applied  to  each  spectral  compon¬ 
ent.  This  analysis  was  carried  out  and  implemented  on  the  computer. 

At  this  point  a  subtle  aspect  of  the  analysis  for  Epstein-type  layers  was  identified. 
All  of  the  profiles  employed  in  all  of  the  available  analytic  solutions  extend  to  infinity. 
Half-space  problems  can  be  treated  by  using  a  profile  that  is  symmetric  about  the 
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origin,  but  the  profile  is  still  infinite  in  one  direction.  In  the  ionospheric  inverse 
scattering  problem,  it  is  necessary  to  have  the  source  a  finite  distance  from  the 
measurement  surface  (e.g.,  the  satellite  must  be  at  finite  altitude).  If  this  is  not  the 
case,  then  only  a  single  plane  wave  illuminates  the  ionosphere.  Furthermore,  the 
scattering  and  measurements  occur  in  the  far  field  of  the  source,  which  implies  that  the 
relative  phase  of  the  unscattered  component  of  the  wave  illuminating  the  ionosphere  at 
various  altitudes  is  the  same.  The  result  is  the  loss  of  the  ability  to  resolve  he  ght 
information  using  a  single  frequency.  Using  the  analytic  profile  method,  the  geometries 
which  can  be  treated  for  a  semi-infinite  profile  with  the  receiver  in  the  free  space 
portion  (such  as  is  the  case  for  the  problem  being  modeled  here)  are  limited  to  those  in 
which  the  source  is  at  infinity,  or  in  which  the  source  is  in  the  same  half-space  as  the 
receiver  (the  backscatter  case).  Truncation  of  the  profile  at  a  finite  distance  renders 
the  analytic  solution  invalid.  It  is  possible  to  obtain  normal  mode  solutions  in  a 
geometry  involving  a  truncated  analytic  profile,  but  the  full,  continuous  spatial 
spectrum  is  required. 

Two  other,  related  techniques  were  also  briefly  considered.  Numerical  ray 
tracing  programs  exist  for  the  ionosphere  with  arbitrary  profiles.  However,  ray  tracing 
is  valid  only  in  the  geometrical  optics  limit.  Solutions  exist  for  propagation  through 
finite  media  consisting  of  discrete  layers.  Unfortunately,  reflections  are  generated 
from  the  artificial  boundaries  between  the  layers.  It  was  recognized  a  little  over  a 
decade  ago  that  this  effect  rendered  such  solutions  for  ionospheric  problems  invalid 
unless  the  layer  thickness  was  small  both  compared  to  the  distance  over  which  the 
profile  varies  and  compared  to  the  wavelength.  This  latter  requirement  makes  the 
number  of  computations  required  impractical  for  the  cases  of  interest  here. 

5.3  SPATIAL  FREQUENCY  DOMAIN  METHODS  FOR  DIRECT  SCATTERING 

SIMULATION 

When  the  assumption,  that  analytic  solutions  for  such  profiles  as  the  Epstein 
profile  could  be  used  to  simulate  the  required  direct  scattering  data,  proved  false,  other 
methods  were  investigated.  The  required  technique  must  be  efficient,  since  substantial 
data  sets  are  involved.  It  must  be  applicable  to  two-  and  three-dimensional  problems, 
and  preferably  to  media  varying  in  two  and  three  dimensions.  It  must  handle  problems 
with  finite  distances  between  source  and  receiver,  and  in  which  forward  scattering  is 
important.  Finally,  techniques  based  on  approximations  to  the  scattering  or  diffraction 
processes  are  less  than  ideal,  since  the  purpose  is  to  test  an  exact  inverse  technique. 
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One  approach  is  to  solve  the  wave  equation  numerically,  using  finite  difference  or 
finite  element  methods.  This  is  probably  not  adequately  efficient,  and  would  be  quite 
difficult  to  implement. 

A  technique  which  satisfies  the  above  criteria  was  developed  by  Bojarski  (Ref  72), 
who  termed  it  the  "k-space"  method.  A  simulation  begins  by  specifying  the  electron 
density  as  a  function  of  the  spatial  coordinates,  and  thus  determining  the  refractive 
index.  It  is  desired  to  compute  the  fields  on  the  ground,  due  to  the  source  of  the 
probing  field,  in  the  presence  of  this  spatially  varying  refractive  index.  To  compute  the 
fields  the  wave  equation,  Equation  2,  must  be  solved.  However,  the  refractive  index 
enters  via  Equation  3,  which  in  turn  depends  on  the  field.  Thus,  these  two  equations 
must  be  solved  simultaneously.  Bojarski  noted  that  in  x  space  Equation  2  is  a 
differential  convolution  equation,  while  Equation  3  is  algebraic.  Conversely,  in  Fourier 
transform  space  (or  "k"  space  in  Bojarski's  notation),  Equation  2  is  algebraic,  while 
Equation  3  is  a  convolution  equation.  He  then  devised  an  iterative  technique  which 
alternates  between  the  two  spaces,  solving  each  equation  in  the  domain  where  it  is 
algebraic.  Because  FFTs  are  used  to  go  back  and  forth  between  the  two  domains,  the 
solution  is  obtained  quite  efficiently.  Often,  only  a  few  (<10)  iterations  are  necessary. 
Kruger  (Ref  73)  applied  this  technique  to  several  two-dimensional  problems. 

The  k-space  technique  was  implemented  by  the  author.  However,  problems  were 
encountered  with  convergence  of  the  iterations  using  realistic  ionospheric  refractive 
index  distributions.  A  simpler,  noniterative  version,  used  previously  by  the  author  to 
simulate  the  Holographic  Radio  Camera  experiment  (Ref  5),  was  then  implemented. 
(After  the  author  developed  it,  it  was  recognized  that  this  method  is  essentially 
equivalent  to  the  first  iteration  of  Bojarski's  technique).  It  models  the  ionosphere  as 
made  up  of  homogeneous  cells,  large  compared  to  a  wavelengtgh.  Propagation  through 
each  cell  is  computed  by  using  the  analytical  spatial  Fourier  transform  of  the  Green's 
function  for  a  homogeneous  medium  with  the  properties  of  the  cell,  multiplying  by  the 
numerical  spatial  Fourier  transform  of  the  field  incident  on  the  cell,  and  computing  the 
inverse  numerical  Fourier  transform  of  the  result.  This  is  a  Fourier  domain  method  of 
evaluating  Equation  10  (the  solution  to  Equation  2)  without  having  to  sample  at  an 
interval  small  compared  to  a  wavelength.  Preliminary  results  indicate  that  this  method 
is  quite  efficient,  and  may  be  accurate  enough  to  provide  the  needed  simulated  data. 

Very  recently  Kastner  and  Mittra  (Refs  74,75)  have  presented  a  modified  version 
of  Bojarski's  k-space  method  which  they  call  the  "spectral  domain"  method.  The 
relationship  between  the  two  methods  can  be  seen  by  considering  a  two-dimensional 
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problem.  In  the  k-space  method,  two-dimensional  Fourier  transforms  are  employed, 
and  consistency  between  Equations  2  and  3  is  enforced  once  per  iteration.  In  the 
spectral  domain  method,  a  series  of  one-dimensional  transforms  is  taken  to  propagate 
the  field  in  the  second  dimension  (as  in  the  author's  approach).  However,  the 
consistency  between  Equations  2  and  3  and  satisfaction  of  ihe  boundary  values  is 
enforced  at  each  cell  boundary — many  times  per  iteration.  The  author  thinks  it  is  likely 
that  this  will  eliminate  the  convergence  problems  mentioned  above.  Discussions  with 
Mittra  (Ref  76)  confirm  that  the  spectral  domain  method  has  been  used  successfully  on 
two-dimensional  inhomogeneous  dielectric  media.  This  technique  should  receive  first 
priority  in  future  work. 

5A  RESULTS  OF  ONE-DIMENSIONAL  TIME  DOMAIN  INVERSE  SCATTERING 

SIMULATION 

As  a  means  of  demonstrating  the  quality  of  electron  density  distribution  recon¬ 
structions  which  should  be  possible  using  inverse  scattering  techniques,  a  one-dimen¬ 
sional  simulation  was  carried  out.  It  was  also  felt  that  the  technique  employed  might 
have  direct  application  to  the  processing  of  data  from  existing  vertical  incidence 
sounders.  This  has  proven  to  be  the  case. 

The  method  used  is  an  efficient  "stepping  in  time"  technique  developed  by 

Bojarski  (Ref  77).  The  equations  are  given  in  his  paper.  As  applied  to  obtain  the  results 

shown  below,  an  impulsive  plane  wave  is  assumed  to  be  incident  vertically  from  below 

on  an  ionosphere  which  has  a  vertical  variation  in  electron  density  such  as  is  shown  in 

Figure  2.  The  horizontal  axis  is  normalized;  the  "F-region"  peak  corresponds  to  a 
8  3 

density  of  3  x  10  electrons/m  .  The  corresponding  value  for  the  bottom  of  the 
"E-region"  ledge  is  7  x  10^  e"/m^,  and  the  value  at  1000  km  altitude  is  3  x  lO^e'VnrA 
This  is  a  reasonably  typical  daytime  profile.  The  time  history  of  the  reflection 
coefficient  is  computed.  Noise  is  added  by  specifying  a  fraction  of  the  peak  magnitude 
of  the  reflection  coefficient,  and  then  for  each  time  sample  computing  a  pseudorandum 
number  in  the  range  of  zero  to  this  fraction  of  the  peak  which  is  added  to  the  reflection 
coefficient  at  that  time.  This  noise-contaminated  reflection  coefficient  is  then 
employed  in  the  inverse  scattering  algorithm  to  reconstruct  the  profile.  The  maximum 
deviation  and  the  root  mean  square  deviation  between  the  true  and  reconstructed 
profile  are  also  computed. 
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Figures  3  through  8  show  the  results  for  various  amounts  of  noise.  Quantitatively, 
these  results  can  be  summarized  as  follows  (single  precision— 6-1/3  digit— computations 
were  used): 


Noise  fraction  of 

Peak  Reflection 
Coefficient 

Maximum 

Deviation 

RMS 

Deviation 

0 

2.98  x 

io-s 

1.08 

x  10"8 

1.0  x  10'2 

1.07  x 

10'2 

4.82 

x  10-3 

5.0  x  10"2 

5.17  x 

10‘2 

2.26 

x  IQ'2 

1.0  x  10_1 

9.75  x 

10'2 

4.17 

x  10‘2 

1.5  x  10’1 

1.36  x 

10"1 

5.75 

x  10-2 

2.0  x  10_I 

1.67  x 

10_1 

7.07 

x  10-2 

Several  things  are  significant  about  these  results.  First,  the  maximum  deviation 
and  the  rms  deviation  are  both  bounded  by  the  noise  fraction.  This  is  consistent  with 
the  results  of  the  author  for  the  Exact  Theory  (see  Section  6).  However,  it  is  in  sharp 
contrast  with  the  results  obtained  and  the  conclusion  drawn,  with  respect  to  stability  in 
the  presence  of  noise,  by  Bojarski  (Ref  77).  The  author  postulates  that  this  is  because 
Bojarski  employed  a  profile  with  a  large  discontinuity  in  it.  It  would  be  valuable  to 
perform  further  numerical  experiments  and  analysis  to  test  this  postulate.  Second,  the 
profile  above  the  maximum  is  reconstructed  from  only  bottom-side  data.  Without 
resorting  to  incoherent  scatter  radar,  the  author  knows  of  no  other  technique  which 
permits  determination  of  the  profile  above  the  peak  using  only  bottom-side  illumination 
and  measurements.  Third,  it  appears  that  this  same  data  processing  can  be  directly 
applied  to  existing  vertical  incidence  sounders.  The  Digisonde  appears  to  be  ideal  for 
this. 


5.5  APPLICABILITY  OF  THE  HOLOGRAPHIC  RADIO  CAMERA  RESULTS  AS  A 
SIMULATION  OF  THE  CURRENT  PROBLEM 

In  Section  6  it  is  shown  that,  in  the  limit  where  the  wavelength  of  the  probing 
radiation  is  short  compared  to  all  distances  in  the  problem  and  compared  to  any  medium 
variations  to  which  the  measurement  is  sensitive,  the  holographic  reconstruction  is 


41 


ALTITUDE  (km) 


RT-21027 


NORMALIZED  ELECTRON  DENSITY 


Figure  3.  Reconstructed  normalized  electron  density  profile  obtained  with  the 
time  domain  inverse  scattering  technique  and  a  noise  fraction  of  zero. 
Maximum  deviation  is  2.98  x  10"*  and  RMS  deviation  is  1.08  x  10"*. 
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NORMALIZED  ELECTRON  DENSITY 
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Figure  6.  Reconstructed  normalized  electron  density  profile  obtained  with  the  time  domain 
inverse  scattering  technique  and  a  noise  fraction  of  1  x  10" l.  Maximum 
deviation  is  9.75  x  10~2  and  RMS  deviation  is  #.17  x  10~2. 
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RT-2I03I  NORMALIZED  ELECTRON  DENSITY 


Figure  7.  Reconstructed  normalized  j _ 

inverse  scattering  technique  and  a  •  ^  Pr°f*le  0bta*ned  wif  **  domain 

deviation  is  1.36  x  10‘J  J^RMS  ^  °‘  ‘  V  10  *  Ma*imum 

1U  *"  RMS  deviation  is  5.75  x  i0~2. 


06 
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Figure  «.  Reconstruct  nonnalit  electron  density  prof  Ue  obtained  -id.  the  tin*  donuun 
inverse  scattering  technique  and  a  noise  fraction  of  2.0  x  10  .  Maxrmum 
deviation  is  1.67  x  10'1  and  RMS  deviation  is  7.07  x  10  . 
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equal  to  ip.  It  is  also  shown  that,  for  reconstructions  of  the  ionosphere,  the  intensity  of 
the  holographic  field  is  proportional  to  the  electron  density.  It  follows  that  the  results 
from  the  author's  Holographic  Radio  Camera  demonstration  experiment,  summarized  in 
Section  6,  are,  in  fact,  excellent  examples  of  the  type  of  reconstruction  obtainable 
using  inverse  scattering  techniques. 


1 

] 
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6.  PRACTICAL  ASPECTS  OF  USING  INVERSE  SCA'  IT-RING  TO 
PROBE  THE  IONOSPHERE 


This  section  considers  several  practical  aspects  of  applying  inverse  scattering 
techniques  to  detemrine  the  three-dimensional  electron  density  distribution  of  the 
ionosphere.  The  effects  of  noise  in  measurements  are  considered  in  Section  6.1. 
Section  6.2  examines  the  effects  of  incomplete,  spatially  sampled  measurements,  and 
the  resolution  which  can  be  achieved.  Since  the  results  in  these  two  subsections  have 
been  available  previously,  they  are  only  summarized  here.  Section  6.3  examines  the 
choice  of  probing  frequency,  and  shows  why  scalar  measurements  may  be  sufficient 
with  the  proper  choice  of  frequency.  Section  6.4  presents  the  important  result  that, 
again  with  the  proper  choice  of  frequency,  and  with  certain  other  conditions  satisfied, 
it  may  not  be  necessary  to  carry  out  the  full  solution  to  the  Exact  Theory  to  obtain  the 
relative  distribution  of  electron  density.  Finally,  Section  6.5  reviews  the  Holographic 
Radio  Camera  demonstration  experiment  as  an  example  of  how  the  type  of  data 
required  for  inverse  scattering  measurements  can  be  taken,  and  what  results  were 
obtained  using  the  simpler  processing  described  in  Section  6.4. 

6.1  THE  EFFECTS  OF  NOISE  MEASUREMENTS 

The  effects  of  measurement  noise  on  the  reconstructed  refractive  index  obtained 
using  the  theory  of  Section  2  can  be  seen  by  writing  p^(x)  =  d^(x),  +  bN  (x)  where  dN(x) 
contains  a  contribution  due  to  noise.  It  follows,  using  Equation  12,  that  this  is 
equivalent  to  a  source  term  p(x)  +  p^Cx),  where  p(x)  is  the  true  source  and  PN(x) 
contains  the  effect  of  the  noise.  From  this  it  can  be  seen  that  the  signal-to-noise  ratio 
of  the  solution  is  the  signal-to-noise  ratio  of  d^(x).  Since  d^(x)  depends  on  the  integral 
over  the  measurement  surface  of  the  measured  field  values  (Equation  4),  the  signal-to- 
noise  ratio  of  P|_j(x)  [and  thus  of  p(x)]  is  not  less  than  (and  may  be  greater  than)  the 
signal-to-noise  ratio  of  the  measured  data.  This  has  been  confirmed  by  numerical 
experiments  (Ref  52).  Recently,  tin  luthor  (Ref  78)  has  proven  analytically  that  the 
signal-to-noise  ratio  of  p(x)  is  greater  than  or  equal  to  the  signal-to-noise  ratio  of  the 
measured  data. 
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6.2  THE  EFFECTS  OF  INCOMPLETE  MEASUREMENTS,  AND  RESOLUTION 

Practical  inverse  scattering  problem  measurements  almost  always  involve  dis¬ 
crete  measurements  over  a  limited  aperture,  as  opposed  to  the  continuous  measure¬ 
ments  over  a  closed  surface  used  in  the  above  theory.  NAager  and  Bleistein  (Ref  79) 
have  shown  that,  in  the  physical  optics  limit,  the  spatial  bandwidth  over  which 
measured  data  is  known  is  the  spatial  bandwidth  over  which  the  source  term  can  be 
determined.  A  similar  but  more  general  result  follows  directly  from  analysis  of  the 
three-dimensional  spatial  Fourier  transform  of  Equations  k  and  12.  Let  be  the  spatial 
frequency  variable,  and  let  capital  letters  denote  the  transformed  functions.  For  the 
general  case,  o^Uf)  is  known  over  the  whole  surface  S,  and  thus  for  all  0_  v_  Vq  (the 
upper  bound  is  rather  than  oo  ,  since  D,  and  thus  S,  are  of  finite  size).  Since  the  left 
side  of  the  transformed  version  of  Equation  12  is  known  for  all  0_  v_  _  v,  it  follows  that 
P(rO  is  determined  over  this  range.  Now  let  <,Mx)  be  measured  at  discrete  points  over  a 
limited  aperture.  Then  O')  is  determined  for  where  these  spatial 

bandlimits  are  determined  by  the  aperture  size  and  sample  spacing.  It  follows  from  the 
transform  of  Equation  4  that  4>^(£)  is  similarly  bandlimited,  and  from  the  transform  of 
Equation  12  that  P(y)  can  be  determined  over  this  band  of  spatial  frequencies. 

A  related  but  more  general  result  on  the  resolution  obtainable  in  a  coherent 

imaging  process  has  been  obtained  by  the  author  (Ref  80).  The  lateral  resolution  for  a 

coherent  process  is  Ax  =  (Az/D)(A<6/27r),  where  X  is  the  wavelength,  z  is  the  distance 

from  the  aperture  to  the  target,  and  D  is  the  size  of  the  aperture  over  which  the  data  is 

recorded.  A  is  the  accuracy  with  which  phase  can  be  measured  across  the  aperture. 

Note  that  as  A d> — ►2w,  the  resolution  goes  to  the  classical  (Rayleigh)  result  for 

incoherent  imaging.  The  second  term,  involvings ,  represents  the  gain  in  resolution 

as  a  result  of  the  added  information  present  in  the  phase  recorded  in  a  coherent 

2 

process.  A  similar  result  is  obtained  for  the  longitudinal  resolution:  A  z  =  (Xz  /D) 
(A<f»/27r).  The  accuracy  with  which  phase  can  be  measured  is  determined  by  the  power 
signal-to-noise  ratio  of  the  measurement,  S/N:  Ad>  =  (2S/N)'  .  These  theoretical 
results  have  been  verified  by  numerical  experiments  (Refs  13,52,80)  and  by  results 
obtained  with  the  Holographic  Radio  Camera  (see  Section  6.5). 
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6.3  THE  CHOICE  OF  PROBING  FREQUENCY:  WHEN  ONLY  SCALAR 

MEASUREMENTS  ARE  NEEDED 

Section  2.4  derived  the  equations  for,  and  discussed  the  solution  to,  the  Exact 
Inverse  Scattering  Theory  for  vector  fields  and  anisotropic  media.  Sections  1.2  and  2.3 
discussed  the  equations  applicable  to  the  description  of  such  m^dia.  The  complications 
introduced  by  anisotropic  media  and  the  need  to  measure  an'1  process  vector  field  data 
are  significant.  It  would  be  very  important  if  conditions  under  which  scalar  field 
measurements  and  scalar  processing  could  yield  electron  density  distributions  could  be 
identified.  Fortunately,  this  is  the  case. 

Consider  Equations  13  and  14  for  the  refractive  index,  n,  of  a  magneto-ionic 

medium.  Assume  that  the  region  of  interest  in  the  ionosphere  is  above  ~  70  km 

altitude.  Now  consider  the  relative  magnitudes  of  the  four  terms  (Equation  14)  which 

contribute  to  n.  Above  ~  70  km,  the  maximum  collision  frequency,  v ,  is  <  5  x  105  Hz 

(c.f.  Ref  2,  Table  1.3).  At  a  probing  frequency  of  >100  MHz,  Z<8x  10_i*.  This  is 

negligible  compared  to  the  1  in  the  denominator  of  Equation  1 3,  and  n  can  be  treated  as 

-5  2 

real.  The  earth's  magnetic  field  has  a  magnitude  of  ~5xl0  web/m  .  Thus 
eB/m<8.8  x  10^,  and  Y,  and  Y-j.  are  <M0~2  at  frequencies  above  100  MHz.  It  follows 
that  terms  of  order  Y2  can  be  neglected  with  respect  to  1  in  the  denominator  of 
Equation  13.  The  result  is  that,  for  frequencies  above  -  100  MHz,  the  refractive  index 
is  given  by  Equation  15.  This  is  independent  of  the  direction  of  propagation  with 
respect  to  the  magnetic  field,  and  the  ionosphere  behaves  isotropically.  At  such 
frequencies,  scalar  measurements  and  scalar  processing  are  sufficient. 

A  frequency  of  100  MHz  or  higher  is  preferred  for  at  least  two  other  reasons.  At 
HF  (3-30  MHz)  frequencies  and  below,  significant  reflection  occurs  in  the  ionosphere. 
Within  the  context  of  inverse  scattering,  this  means  that  the  backscattered  component 
becomes  significant,  and  reconstructions  using  measurements  primarily  sensitive  to  the 
forward  scattered  component  (e.g.,  ground-based  measurements  with  a  satellite-borne 
source)  would  suffer  a  loss  of  sensitivity.  Also,  structure  of  interest  may  have  scale 
sizes  of  the  order  of  from  several  tens  of  meters  to  100  meters.  As  discussed  in 
Section  6.4,  it  is  often  desirable  to  have  the  wavelength  be  small  compared  to  the 
structure. 

6.4  THE  RELATIONSHIP  BETWEEN  HOLOGRAPHY  AND  INVERSE  SCATTERING 

It  was  noted  in  Section  2.1  that  the  quantity  in  which  the  field  measurements 
appear,  <f>^  (Equation  4),  is  the  field  obtained  by  reconstructing  a  hologram  recorded 
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over  the  measurement  surface.  It  is  therefore  to  be  expected  that  there  would  be  a 
close  relationship  between  holography  and  inverse  scattering.  Indeed,  this  is  the  case, 
und  a  very  important  approximation  comes  out  of  this  relationship. 

Optical  holographers  often  seek  to  probe  an  inhomogeneous  medium  by  obtaining 
images  at  various  depths.  Indeed,  humans  solve  an  inverse  scattering  problem  every 
time  they  recognize  an  object  from  its  image.  Imaging  is  an  inverse  scattering 
problem:  Based  on  observed  field  data,  information  about  the  scatterer  (contained  in 
p(x)  in  the  formulation  of  Section  2.l)  is  sought.  If  only  the  shape  of  the  scatterer  is 
sought,  then  it  can  be  shown  that  a  physical  optics  far-field  approximation  to  the  exact 
inverse  scattering  solution  is  often  sufficient  (Refs  11,81).  However,  if  any  information 
other  than  the  shape  is  desired  (e.g.,  the  complex  refractive  index),  then  the  full  exact 
solution  is  required. 

It  was  noted  above  that<£^  is  the  field  obtained  from  a  holographic  reconstruc¬ 
tion.  From  Equation  9  , 

4>H  =  <»  -y*dv  g *P  (51) 

Several  conclusions  follow  from  this  equation: 

1.  <f>  is  a  field.  It  is  not  the  source  term,  p.  The  source  term,  p,  contains  all 

the  information  about  the  source.  Simple  examination  of  the  holographic 

reconstruction,  cannot  in  general  give  p .  To  obtain  p  ,  Equation  12  must 
be  solved. 

2.  Holographers  commonly  make  the  statement,  "The  holographic  reconstruction 
is  the  field  which  existed  in  object  ^sourcej  space."  This  is  an  approximation 
which,  from  Equation  51,  is  given  by 

4>H  =  (52) 

3.  If,  and  only  if,  the  scatterer  is  of  such  a  form,  and  if,  and  only  if,  the 
wavelength  and  all  other  dimensions  in  the  problem  are  such  that  the  field 
scattered  by  the  scatterer  has  a  perturbation  introduced  into  it  which  is  the 
shape  of  the  scatterer,  then,  to  the  degree  that  Equation  52  holds,  it  is 
possible  to  deduce  the  shape  of  the  scatterer  from  a  holographic 
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reconstruction.  In  other  words,  the  shape  can  be  deduced  when  the  field 
"looks  like"  the  scatterer.  This  is  often  true  for  objects  with  sharp  edges 
imaged  using  optical  holography. 

4.  The  approximation  of  Equation  52  can  only  be  true  when  the  second  term  of 
Equation  51  is  0.  The  kernel,  g*,  is  of  the  form  |cos(kr)  -  i  sin(kr)j/47rr.  The 
integral  will  thus  be  small  as  r-»oo.  Thus,  the  statements  of  paragraphs  ?  and 
3  hold  only  in  the  far  field.  Note  that  this  requires  that  reconstructions, 
"near"  enough  the  position  of  the  scatter er  so  that  the  feld  "looks  like"  the 
scatterer,  must  still  be  in  the  far  field  in  terms  of  wavelengths. 

The  statement,  that  the  requirement  that  the  scattered  field  "look  like"  the 
scatterers,  is  often  true  for  objects  with  edges  that  are  "sharp"  and  can  be  expressed  in 
a  more  quantitative  manner.  Consider  the  kernel  of  the  integral  in  Equation  12.  In 
three  dimensions,  Img  is  given  by 


Img  =  sin(kr)/4-rrr 

(53) 

Using  convolution  notation,  Equation  12  can  be  written  as 

</>  j_j  =  2i  p  *  sin(kr)/477-r 

(54) 

Now  consider  the  short  wavelength  limit.  As  X— *0,  sin(kr)/4« 
a  three-dimensional  delta  distribution),  and 

r  ->5(r)/2  (where  6  denotes 

lim  <J>H  =  2ip*6(r)/2 

X— 0 

(55) 

lim  <f>H  =  ip 

(56) 

X  — *0 

Physically,  this  says  that  when  the  wavelength  is  small  enough  compared  to  any 
structure  in  the  scatterer  (i.e.,  the  edges  of  the  scatterer  are  "sharp"),  the  convolution 
with  Img  acts  like  a  convolution  with  a  delta  distribution,  replicating  the  scatterer  in 
the  field  <£H. 

Equation  56  is  a  very  important  result  for  ionospheric  probing.  Consider  the  use 
of  frequencies  above  100  MHz,  and  restrict  interest  to  spatial  resolution  large  com¬ 
pared  to  the  wavelength.  Then  Equation  56  is  a  good  approximation,  and  Equation  15 
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for  the  refractive  index  also  applies.  From  Equation  3  it  can  be  seen  that  p  (x)  is 
proportional  to  n  (x)-l,  which  in  turn,  from  Equa'  •  15,  is  proportional  to  N(x),  the 
desired  spatial  distribution  of  electron  de  , .  Suppose  that  the  quantity 
fcf>^(x)<i>H  (x)  --the  magnitude  of  the  reconstructed  holographic  field — is  computed. 
From  Equation  56  this  will  be  proportional  top  (x)  and  thus  to  N(x).  Providing  the  short 
wavelength  limit  is  satisfied,  the  magnitude  of  the  reconstructed  holographic  field  is 
proportional  to  the  electron  density  distribution. 

One  of  the  major  reasons  this  is  such  a  powerful  result  is  the  accuracy  and 
efficiency  with  which  4>^(x)  can  be  computed,  even  in  the  presence  of  noise.  Let  f  and 
denote  spatial  Fourier  direct  and  inverse  transforms,  respectively,  computed  using 
FFT  techniques.  Then  it  can  be  easily  shown  (Ref  82,  Sections  3-7;  Ref  83,  Section  5.4 
and  Appendix  I)  that: 


0H(x,y,z) 


=  r 


j/  (x,y,0)|exp|ikz|^l-(\yx)2-( 


(57) 


In  Equation  57,  </>  (x,y,Q)  is  the  data  recorded  on  the  ground,  v  and  v  are  spatial 

x  y 

frequency  variables,  and  the  Fourier  transform  operations  are  two  dimensional,  from 

(x,y)  space  to  ( v  ,  v J  space.  Equation  57  is  derived  from  Equation  4  by  noting  that  the 
x  y 

integral  in  Equation  4  involves  two-dimensional  spatial  convolutions,  applying  the 
convolution  theorem  of  Fourier  analysis,  and  analytically  evaluating  the  two- 
dimensional  spatial  Fourier  transform  of  g*.  Direct  evaluation  of  Equation  4  using 
numerical  integration  would  require  sampling  the  integrand  at  intervals  small  compared 
to  the  periods  of  its  oscillation.  These  periods  are  of  the  order  of  \/d,  where  d  is  any 
dimension  of  interest.  This  is  clearly  impossible  from  either  a  measurement  or  a 
computation  standpoint.  Equation  57  totally  avoids  this  problem,  and  is  more  efficient 
and  accurate  as  well. 


6.5  THE  HOLOGRAPIC  RADIO  CAMERA 

The  previous  subsection  explained  the  close  relationship  between  holography  and 
inverse  scattering.  The  Holographic  Radio  Camera  demonstration  experiment 
(Refs  4,5)  was  carried  out  to  provide  a  visualization  of  ionospheric  structure:  Images  of 
the  electronic  density  distribution.  A  summary  of  the  experiment,  its  results,  and  the 
important  insight  about  ionospheric  structure  obtained  from  these  results  is  given  in 
this  subsection.  This  material  is  included  in  this  paper  because  it  is  an  example  of  how 
the  data  used  in  inverse  scattering  can  be  recorded,  and  an  excellent  demonstration  of 
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the  practicality  of  such  measurements.  In  particular,  it  demonstrates  how  images  of 
the  ionospheric  electron  density  distribution  can  be  obtained  using  the  simplified 
inverse  scattering  techniques  of  the  previous  subsection. 

Historically,  the  ionosphere  has  been  viewed  as  having  a  structure  dominated  by 
random  turbulence.  The  primary  basis  for  this  picture  r.as  been  ground-based 
observations  of  the  temporal  spectra  of  the  fluctuations,  or  scintillations,  introduced  on 
signals  from  satellites  or  stellar  radio  sources.  Because  of  the  lateral  motion  of  the 
ionosphere  with  respect  to  the  source,  these  temporal  spectra  can  be  interpreted  in 
terms  of  the  spatial  frequency  spectrum  of  the  irregularity  structure  in  the  medium. 
These  spectra  display  a  power-law  dependence  of  power  spectral  density  on  spatial 
frequency.  Because  this  is  consistent  with  the  spectra  which  would  be  observed  if  the 
ionosphere  had  a  turbulence-dominated  structure,  such  a  structure  has  been  inferred. 
The  results  thus  far  from  the  Holographic  Radio  Camera  suggest  that  a  far  more 
orderly  structure  may,  in  fact,  be  the  case,  at  least  for  the  mid-latitude  ionosphere. 

Holography  is  the  process  of  recording,  coherently,  *he  amplitude  and  phase  of  a 
wave.  The  Holographic  Radio  Camera  accomplishes  this  at  a 'm  wavelength.  In  the 
demonstration  experiment  described  here,  32  antennas,  separated  32  m  apart,  were  laid 
out  along  an  east-west  line  (near  La  Posta,  California,  at  a  latitude  of  -~32°N).  The 
same  local  oscillator  signal  was  distributed  coherently  along  this  1  km  line  array,  and 
the  intermediate  frequency  signals  were  sampled,  digitized,  and  recorded  such  that  the 
phase  and  amplitude  of  the  field  at  each  antenna  was  recorded  coherently  as  a  function 
of  time.  The  radio  sources  for  the  experiment  were  the  Navy  Navigation  Satellite 
System  (NNSS)  satellites,  which  are  in  circular  polar  orbits  at  an  altitude  of  1000  km. 
For  orbits  on  which  the  satellite  was  at  a  high  enough  elevation  angle  to  be  in  the  main 
lobe  of  the  vertically  pointing  antennas,  the  motion  of  the  satellite,  during  the  portion 
of  the  pass  when  it  was  directly  overhead,  was  approximately  perpendicular  to  the  line 
of  the  array.  By  sampling  in  time,  aperture  synthesis  occurred  and  a  series  of  ~  1  km 
by  ~1  km  radio  frequency  holograms,  each  consisting  of  32  x  32  points,  was  obtained. 
The  details  of  the  sampling  and  recording,  and  an  in-depth  analysis  of  the  aperture 
synthesis  approximation  and  of  the  resulting  corrections  applied  to  the  data  are  given  in 
Reference  4. 

The  resulting  holograms  were  processed  by  a  computer  to  obtain  reconstructions 
of  the  ionosphere  at  a  series  of  altitudes  between  200  km  and  600  km.  Mathematically, 
the  reconstruction  process  involved  evaluating  Equation  57.  Because  of  the  choice  of 


wavelength  and  the  resolution  and  the  geometry  of  the  measurements,  the  field 
magnitudes  presented  below  can  be  considered  to  be  proportional  to  the  electdron 
density  distribution  at  the  altitude  of  reconstruction  (see  Section  6 .4). 

6.5.1  Experimental  Results 

Figures  9  through  12  show  examples  of  the  reconstructions  obtained  to  date  with 
the  Radio  Camera.  The  data  for  these  figures  were  taken  on  Ouly  24,  1976,  at 
approximately  5:19  a.m.  local  time.  The  satellite  passed  over  the  array  within  three 
degrees  of  the  north-south  line.  The  two  end  antennas  were  recorded  as  having  zero 
signal  due  to  cable  problems.  The  data  frames  used  to  produce  these  figures  were  taken 
over  approximately  a  30  second  time  span  centered  about  the  zenith  time  (5:19  a.m.). 
Vertical  incidence  'onograms  recorded  adjacent  to  the  array,  five  minutes  before  and 
after  the  zenith  time,  indicate  that  the  ionosphere  was  "quiet"  and  relatively  unchang¬ 
ing,  with  some  slight  sporadic  E  present. 

The  field  magnitude  reconstructed  at  an  altitude  of  300  km  is  shown  in  Figure  9. 
The  east-west  and  north-south  directions  are  as  indicated  and  this  orientation  of  the 
reconstructions  is  used  in  all  the  figures.  The  sampling  interval  is  32  m  east-west,  and 
approximately  24.5  m  north-south.  Figure  10  shows  a  number  of  reconstructions  of  the 
field  magnitude  from  a  single  frame  of  data,  at  the  indicated  altitudes.  The  0  km 
diagram  represents  the  magnitude  of  the  field  recorded  on  the  ground  by  Radio  Camera. 
The  striking  feature  present  in  both  figures  is  the  highly  periodic  structure,  present  in 
both  directions.  At  300  km  altitude,  the  period  of  the  structure  is  about  800  m  north- 
south  and  about  300  m  east-west.  As  shown  in  Figure  10,  the  periods  in  both  directions 
remain  approximately  constant  over  the  range  of  altitudes  from  300  to  5 00  m.  There  is 
significantly  more  fine-period  structure  present  at  the  200  m  altitude  than  for  higher 
altitudes. 

Figure  1 1  shows  the  reconstructions  of  the  phase  of  the  field  which  correspond  to 
the  magnitude  data  in  Figure  10.  The  dimensions  are  the  same  for  both  sets  of 
reconstructions.  The  "vertical"  axis  on  each  reconstruction  is  the  same,  representing  a 
maximum  phase  excursion  from  -  it  to  tt  radians.  In  order  to  permit  the  structure  to  be 
seen  from  a  single  viewing  direction  in  the  isometric,  hidden  line  display,  noise  peaks  in 
the  phase  reconstructions  have  been  zeroed  (this  is  done  for  display  purposes  only). 
Note  that  the  same  type  of  periodocity  present  in  the  magnitude  is  seen  in  the  phase. 
However,  there  is  a  noticeable  decrease  in  the  magnitude  of  the  phase  excursions  near 
the  altitude  of  400  km. 
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Figure  10.  Reconstructions  of  the  magnitude  of  the  field  at  the  indicated  altitudes 
The  effective  dimensions  of  the  array  in  each  direction  are  shown  for 
each  altitude. 


Figure  1 1.  Reconstructions  of  the  phase 
of  the  field  at  the  indicated 
altitudes 


Figure  12.  Reconstructions  ol  the  field  magnitude  at  the  indicated  altitudes  for  successive 
frames  of  data.  The  spacing  between  frames  at  a  given  altitude  is  approxi¬ 
mately  6^5  times  the  north-south  dimension  of  the  reconstruction  at  that 
altitude,  as  shown  in  Figure  9.  The  data  for  these  reconstructions  were 
recorded  on  3uly  2b,  1976  at  3(19  ajn.  local  time. 
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The  method  of  visually  displaying  the  data  actually  used  in  the  laboratory  permits 
the  reconstruction  to  be  oriented  at  an  arbitrary  angle  in  three-space  with  respect  to 
the  viewer,  and  with  either  isometric  or  full  perspective  display.  Considerable 
information  is  lost  by  fixing  the  orientation  and  reducing  the  plots  for  publication.  As  a 
result,  some  of  the  geometrical  aspects  mentioned  in  this  sectijn  may  appear  to  be  less 
convincing  in  the  printed  figures  than  they  are  in  the  laborato  y. 

Figure  12  presents  magnitude  reconstructions  from  twelve  frames  of  data.  The 
satellite  was  at  zenith  between  the  sixth  and  seventh  frames  from  the  left.  Note  that 
there  is  a  significant  amount  of  uniformity  in  the  periodicity  of  the  structure  at  a  given 
altitude  in  each  direction.  Furthermore,  the  period  of  the  structure  in  each  direction  is 
approximately  the  same  at  300,  400,  and  500  km  altitude  for  many  frames.  Much  finer 
structure  is  observed  at  200  km  in  both  directions. 

Observations  on  several  other  days,  at  different  local  times,  produced  similar 
(although  by  no  means  identical)  results. 

6.5.2  Interpretation 

For  the  ionosphere,  significant  levels  of  electron  density  often  exist  at  all 
altitudes  between  ~80  km  and  the  altitude  of  the  satellite.  However,  it  is  precisely  an 
assumption  of  free  space  propagation  upon  which  the  concept  of  modeling  the 
ionosphere  as  a  "thin  screen"  is  based  ^Booker,  Ratcliffe,  and  Shinn  (Ref  84);  Rufenach 
(Ref  85);  Bramley  (Ref  86)J.  In  this  model  the  fluctuations  imposed  on  the  radio  wave 
due  to  passage  through  the  ionosphere  are  confined  to  a  thin  layer  (typically  assigned  an 
altitude  in  the  range  from  200  to  400  km).  Immediately  below  this  layer  the 
fluctuations  are  only  in  the  phase  of  the  wave.  As  the  wave  propagates  through  free 
space  to  the  ground,  amplitude  fluctuations  develop.  Thus,  the  reconstructions  of  the 
field  at  a  given  altitude  as  presented  above  represent  the  field  which  would  have  to  be 
postulated  to  exist  just  below  a  screen  at  that  altitude  to  produce  the  observed  field  on 
the  ground.  Note  that  in  order  for  the  screen  to  be  a  phase  screen,  the  reconstructed 
field  should  have  essentially  uniform  amplitude:  AA/A  should  be  small  compared  to 
A  </>  (where  A  is  the  amplitude  and  $  is  the  phase).  If  this  were  the  case,  the  plot  of  the 
magnitude  of  the  field  would  appear  as  a  smooth  plane  at  the  altitude  of  the  screen, 
while  the  phase  plot  would  show  significant  fluctuations.  None  of  the  frames  of  data 
shown  above  or  obtained  with  the  Radio  Camera  display  this  property  at  the  altitudes 
shown,  between  200  and  500  km.  The  conclusion  is  that  this  experimental  data  cannot 
be  explained  in  terms  of  a  thin  phase  screen  within  this  altitude  range. 
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The  reconstructions  obtained  to  date  have  a  single,  visually  striking  feature  in 
common — the  predominate  periodic  structure.  Although  the  currently  available  sample 
of  processed  reconstructions  is  too  small  to  support  firm  conclusions  about  the 
mechanism  producing  this  structure,  it  is  interesting  to  speculate  on  the  possibility  that 
waves  in  the  ionospheric  medium  could  be  this  mechanism.  The  primary  setni- 
quantitative  observable  involved  is  the  apparent  wavelength.  However,  the  spaces 
between  frames  of  data  in  the  north-south  direction  and  the  resolution  limiations 
prevent  accurate  determination  of  the  wavelength,  and  make  determination  of  the 
horizontal  distance  over  which  the  periodic  structure  extends  uncertain.  Until 
additional  data  is  available,  interpretation  of  the  results  will  be  limited  to  consideration 
of  whether  structure  such  as  that  present  in  the  reconstructions  is  consistent  with  the 
results  of  other  experiments. 

Current  practice  is  to  model  the  ionospheric  plasma  as  a  largely  random, 
turbulent  medium.  The  experiments  to  date  upon  which  the  most  common  pictures  of 
ionospheric  structure  are  based  use,  almost  without  exception,  measurements  inte¬ 
grated  along  the  line  of  sight  from  the  satellite  to  the  ground  (e.g.,  scintillation 
measurements).  The  most  common  analysis  techniques  applied  to  such  measurements 
involve  studies  of  the  signal  statistics,  making  an  instantaneous  description  of  the 
structure  underlying  a  measured  signal  difficult.  Spaced  receiver  measurements  and 
analyses  of  the  spectra  of  scintillations  do  permit  the  irregularity  size  to  be  deter¬ 
mined.  The  dominant  range  of  irregularity  sizes  reported  for  such  measurements,  at 
mid-latitudes,  is  from  MOO  m  to  l  to  2  km  (this  upper  bound  is  probably  due  only  to 
measurement  limitations).  These  scales  are  of  the  order  of  the  spatial  periods  present 
in  the  Radio  Camera  reconstructions. 

However,  it  is  not  immediately  obvious  that  these  Radio  Camera  results  are 
consistent  with  the  fluctuation  power  spectra  typically  recorded  in  VHF  scintillation 
experiments  in  mid-latitudes.  The  power  spectra  from  individual  frames  of  the  Radio 
Camera  data  show  peaks  corresponding  to  the  spatial  periods  present  in  the  reconstruc¬ 
tions.  Typical  results  from  the  VHF  scintillation  measurements  of  other  authors  involve 
measurements  over  periods  of  time  long  compared  to  the  time  in  which  a  frame  of 
Radio  Camera  data  was  recorded.  These  results  from  other  experiments  show  an 
inverse  power  law  spectrum  of  irregularities  which  has  been  used  to  infer  the  turbulent 
medium  model.  When  spectra  from  several  successive  frames  of  Radio  Camera 
observations  are  averaged,  similar  inverse  power  law  spectra  are  also  obtained.  This 
suggests  that  the  difference  between  the  two  types  of  spectra  might,  in  large  measure, 
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be  due  to  differences  in  the  experiments  in  which  they  were  recorded.  The  Radio 
Camera  provides  an  instantaneous  "snapshot"  of  the  diffraction  pattern,  whereas  the 
other  techniques  involve  a  lot  of  averaging.  It  seems  possible  that  the  Radio  Camera  is 
revealing  a  process  of  nonlinear  interaction  between  different  Fourier  components  of 
the  plasma,  which  interaction  in  the  long  run  leads  statistically  to  an  inverse  power  law 
spectrum  of  plasma  turbulence.  In  other  words,  what  has  heretofore  been  described 
only  as  turbulence  may,  in  fact,  have  an  underlying  wave-like  structure.  When  the 
effects  of  propagation  through  this  structure  are  averaged  both  spatially  and  tempor¬ 
ally,  the  structure  appears  turbulent.  Indeed,  Booker  (Refs  87,88)  has  put  forth  a 
consistent  model  of  a  whole  range  of  ionospheric-radio  wave  interactions  in  which  the 
underlying  mechanism  is  nonlinear  coupling  of  acoustic  gravity  wave  energy  (generated 
in  the  neutral  atmosphere)  into  a  spectrum  of  plasma  waves. 


7.  CONCLUSIONS  AND  RECOMMENDATIONS 


7.1  CONCLUSIONS 

The  following  conclusions  are  drawn  from  this  research.  The  conclusions  listed 
are  those  that  are  new,  based  on  the  research  reported  in  this  paper.  Conclusions 
available  from  review  material  included  in  this  paper  are  not  listed.  These  conclusions 
are  listed  in  their  order  of  presentation  in  previous  sections. 

1.  The  Exact  Inverse  Scattering  Theory  is  the  best  available  technique  for 
determining  the  three-dimensional  electron  density  distribution  of  the 
ionosphere,  using  measurements  of  the  scattering  of  electromagnetic  waves 
entering  the  ionosphere  from  outside  the  earth's  atmosphere. 

2.  The  Exact  Inverse  Scattering  Theory  for  vector  fields  has  oeen  presented.  It 
has  been  shown  how,  within  the  context  of  the  magneto-ionic  theory, 
measured  vector  field  data  can  be  used  to  obtain  the  three-dimensional 
distribution  of  both  electron  density  and  collison  frequency  in  the  presence  of 
the  earth's  magnetic  field. 

3.  A  method  of  processing  vector  field  data  using  the  scalar  equations  of  the 
Exact  Theory  to  obtain  the  electron  density  and  collision  frequency 
distributions  was  obtained. 

4.  The  uniqueness  of  the  solution  to  the  inverse  scattering  problem  has  been 
proven  for  the  vector  electromagnetic  field  case. 

5.  Three  additional  proofs  have  been  presented,  substantiating  the  previous 
proof  by  the  author  that  nonradiating  sources  of  compact  support  are 
physically  inconsistent  with  Maxwell's  equations,  and  thus  do  not  exist. 

6.  An  expression  for  the  fundamental  quantity,  in  Bojarski's  exact,  closed-form, 
analytic  solution  for  the  total  field  in  the  inverse  scattering  problem,  in 
terms  only  of  Fourier  transform  operations  on  the  measured  fields  was 
presented. 


7.  An  exact,  closed-form,  analytic  solution  for  the  source  term  in  the  inverse 
scattering  problem  has  been  obtained. 

8.  The  numerical  solution  to  the  exact  inverse  scattering  integral  equation  was 
implemented  and  studied  for  the  two-dimensional  c  lse.  Several  preliminary 
results  related  to  this  numerical  solution  and  to  deconvolution  problems  in 
general  were  obtained,  as  given  in  Section  5.1. 

9.  It  was  shown  that  analytic  solutions  of  the  type  associated  with  many 
solutions  to  the  wave  equation  in  terms  of  special  functions,  for  analytic 
refractive  index  profiles  (e.g.,  the  Epstein  profile),  are  not  applicable  to 
modeling  forward  scattering  ionospheric  experiments  of  the  type  considered 
in  this  paper.  This  is  a  consequence  of  the  inherent  infinite  or  semi-infinite 
extent  of  such  profiles. 

10.  Several  spatial  frequency  domain  techniques  for  simulating  the  problem 
considered  in  this  work  were  implemented,  although  testing  was  not 
completed.  It  is  felt  that  a  useful  approach  was  identified. 

11.  A  one-dimensional,  time  domain  inverse  scattering  technique  was  simulated, 
including  the  effects  of  noisy  data.  The  technique  was  demonstrated  to  be 
stable  with  respect  to  noise  for  typical  ionspheric  profiles.  Errors  in  the 
reconstructed  profile  were  less  than  or  equal  to  the  noise  in  the  data, 
contrary  to  previous  results.  The  technique  permits  the  profile  at  altitudes 
above  the  maximum  to  be  obtained  using  only  bottomside  vertical  incidence 
sounder  data.  The  approach  appears  readily  applicable  to  the  processing  of 
data  from  existing  sounders. 

12.  It  has  been  shown  that  by  choosing  a  short  enough  probing  wavelength,  the 
anisotropic  properites  of  the  ionosphere  are  negligable  within  the  magneto¬ 
ionic  theory  approximations.  This  permits  the  use  of  scalar  field 
measurements  and  scalar  processing  to  obtain  an  accurate  reconstruction  of 
the  electron  density  distribution. 

13.  It  has  been  shown  that  in  the  short  wavelength  limit,  the  field  reconstructed 
from  a  hologram  is  equal  to  i  times  the  source  term  in  the  inverse  scattering 
problem.  Furthermore,  for  the  case  of  the  ionosphere  within  the  framework 


of  the  magneto-ionic  theory,  it  has  been  shown  that  the  magnitude  of  the 
field  as  a  function  of  x  obtained  from  a  holographic  reconstruction  is  directly 
proportional  to  the  electron  density  as  a  function  of  x. 

7.2  RECOMMENDATIONS 

The  above  conclusions  each  provide  numerous  possibilities  for  further  research. 
The  following  recommendations  highlight  those  areas  deemed  to  be  of  most  importance 
in  the  application  of  inverse  scattering  to  the  determination  of  the  ionospheric  electron 
density  distribution. 

1.  Numerical  evaluation  of  the  exact,  closed-form,  analytic  solution  for  the 
source  should  be  studied,  and  the  properties  of  the  solution  examined. 

2.  The  deconvolution  technique  identified  as  most  promising  in  Section  5.1 
should  be  investigated,  both  analytically  and  numerically. 

3.  The  simulation  of  the  determination  of  the  ionospheric  density  distribution 
should  be  developed,  utilizing  the  spatial  frequency  domain  techniques 
identified  in  this  research  to  model  the  direct  scattering. 

4.  The  one-dimensional,  time  domain  inverse  technique  should  be  studied, 
especially  with  regard  to  its  properties  in  the  presence  of  noise  in  the 
measurements.  Application  of  the  technique  to  data  from  existing  sounders 
should  be  pursued. 

5.  An  experimental  observatory,  using  measurements  such  as  were  made  in  the 
Holographic  Radio  Camera  demonstration  experiment,  processed  using  the 
inverse  scattering  techniques  of  this  paper,  should  be  established. 
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